AD  No. AbAb55898 

OOC  FILE  COPY 


Report  No;  FAA-RD -77-189 J 


GENERAL  AVIATION  AIRPLANE 
STRUCTURAL  CRASHWORTHINESS  USER'S  MANUAL 


VOLUME  I 

PROGRAM  "KRASH"  THEORY 


Max  A.  Gamon 


February  1978 


Document  is  available  to  the  U.S.  public  through 
the  National  Technical  Information  Service, 
Springfield,  Virginia  22161. 


Prepared  for  . 

U.S.  DEPARTMENT  OF  TRANSPORTATION 


FEDERAL  AVIATION  ADMINISTRATION 
Systems  Research  A Devetopment  Service 
Washington,  D.C.  20590 


78  0 7 08 


1 


J 


NOTICE 


This  document  is  disseminated  under  the  sponsorship  of 
the  Department  of  Transportation  in  the  interest  of  information 
exchange.  The  United  States  government  assumes  no  liability 
for  its  contents  or  use  thereof. 


DISCLAIMER  NOTICE 


THIS  DOCUMENT  IS  BEST  QUALITY 
PRACTICABLE.  THE  COPY  FURNISHED 
TO  DDC  CONTAINED  A SIGNIFICANT 
NUMBER  OF  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


is  doctiment  provides  a comprehensive  description  of  program  KRASH,  as 
modified* ‘under  Contract  DOT-FA75-WA-3707.^ Included  in  this  User's  Manual  are  the^ 
following  sections. 

; _ -4 

1. '^  Program  KRASH^{pescription  and  Theorj^^'  ' 

2.  ^Program  Controls^  and  Listing,-. 

^■J^e  User's  Manual  is  part  one  of  a three-volume  docximent.  The  second  volume 
describes  modeling  techniques  applicable  to  program  KRASH  to  assist  in  the  develop- 
ment of  improved  structural  crashworthiness  designs , as  well  as  the  detailed  input 
requirements  and  output  capabilities  of  KRASH.  Volume  III  contains  general  back- 
ground data  regarding  structural  crashworthiness  design Each  volume  has  been 
established  in  such  a manner  that  they  can  readily  be  upqated  as  more  data  becomes 
available.  The  subject  matter  contained  within  each  section  can  be  expanded  or 
revised,  as  necessary,  without  affecting  the  other  sections.  Each  section  contains 
its  own  numbering  system  which  facilitates  the  task  of  updating  the  document. 


\ 

t7.  K«y  Words 

16*  Distribution  Stotsmont 

Program  KRASH,  Theory,  Program  Controls, 

Docimient  is  available  to  the  public 

Listing,  Math  Model,  Forces,  Deflections, 

through  the  National  Technical 

Energy,  Stress,  DRI,  Accelerations, 

Information  Service,  Springfield,  VA  22151 

Moments,  Velocities,  Symmetrical, 

Unsyramet rical 

19.  Soourify  Clossif.  (o^  this  roport) 

20.  Socurity  Ctossif.  (of  this  pofo) 

21*  No.  of  Pogo% 

22.  Prie. 

Unclassified 

Unclassl fled 

272 

fin  DOT  F 1700.7  (8-72) 


It  •production 


4 on 


U i 


FOREWORD 


This  report  was  prepared  by  the  Lockheed-California  Company  xinder 
contract  DOT-FA75-WA-370T.  The  report  contains  a partial  description  of  the 
effort  performed  as  Task  II  of  a three  task  effort.  The  Task  II  report  covers 
the  period  from  July  1976  to  December  1977.  The  work  was  administrated  xinder 
the  direction  of  the  Federal  Aviation  Administration  with  H.  Spicer  acting  as 
Technical  Monitor. 

The  project  leader  was  Gil  Wittlin  of  the  Lockheed-CeJ-ifornia  Company. 
Important  contributions  were  made  to  the  program  by  the  Cessna  Aircraft  Com- 
pany under  the  direction  of  D.J.  Ahrens  and  W.B.  Bloedel.  M.A.  Gamon  and 
W.  L.  LaBarge  participated  in  the  development  of  the  User's  Manual.  The 
Lockheed  effort  was  supervised  by  J.E,  Wignot  (Dynamics  Loads  Group)  and 
R.F.  O'Connell  (Aeromechanics  Department). 
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This  docvunent  provides  a con^jrehensive  description  of  program  KRASH,  as 
modified  under  Contract  DOT-FA75-WA-3707.  Included  in  this  User's  Manual 
I are  the  following  sections: 


1.  Program  KRASH;  Description  smd  Theory 

2.  Progrsim  Controls  and  Listing 


The  User's  Manual  is  part  one  of  a three-volume  document.  The  second 
volume  describes  modeling  techniques  applicable  to  program  KRASH  to  assist 
in  the  development  of  improved  structural  crashworthiness  designs,  as  well 
as  the  detailed  input  requirements  and  output  capabilities  of  KRASH.  Voliune 
III  contains  general  background  data  regarding  structiiral  crashworthiness 
design.  Each  volume  has  been  established  in  such  a manner  that  they  can 
readily  be  updated  as  more  data  becomes  available.  The  subject  matter  con- 
tained within  each  section  can  be  expanded  or  revised,  as  necessary,  with- 
out affecting  the  other  sections.  Each  section  contains  its  own  numbering 
system  which  facilitates  the  task  of  updating  the  document . 
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SECTION  1 


PROGRAM  KRASH 


1.1  INTRODUCTION 

This  section  presents  a comprehensive  description  of  the  crash  analysis  digi- 
tal computer  program  called  KRASH.  Section  1.2  gives  an  outline  of  the  general 
features  and  capabilities  of  the  KRASH  program.  Section  1.3  presents  in  detail 
all  the  theoretical  development  of  the  equations  contained  in  KRASH.  Section  2 
outlines  the  input-output  control  parameters  in  KRASH,  and  presents  a siir^jlified 
program  flow  diagram.  A coniplete  program  listing  is  included  in  Appendix  A. 

1.2  PROGRAM  KRASH  DESCRIPTION 

Digital  computer  program  KRASH  predicts  the  response  of  vehicles  to  multi- 
directional crash  environments.  The  program  computes  the  time  histories  of  N 
interconnected  masses,  each  allowed  six  degrees  of  freedom  defined  by  inertial 
coordinates  x^^,  y^ , and  Eulerian  angles  <j)^,  0^,  where  i = 1,  2...N. 

Euler's  equations  of  motion  are  written  for  each  mass.  The  equations  of  motion 
are  Integrated  numerically  to  obtain  velocities,  displacements,  and  rotations. 
Gravitational  forces,  internal  forces  and  moments,  and  external  forces  are 
computed.  For  small  deflections  a linear  analysis  is  obtained,  and  for  large 
deflections,  general  plastic  deformation  is  allowed.  The  program  provides  for 
unloading  and  subsequent  reloading  along  a linear  elastic  line. 

Program  KRASH  describes  the  interaction  between  a series  of  massless  inter- 
connecting structural  elements  and  concentrated  rigid  body  masses  to  which 
the  structural  elements  are  attached  at  their  ends  with  the  appropriate  end 
fixity  (pinned,  fixed).  The  structural  elements  can  be  connected  between 
"node  points"  which  are  offset  from  and  rigidly  attached  to  selected  mass 
points.  The  interconnecting  elements  represent  the  stiffness  characteristics 
of  the  structure  between  the  masses.  The  masses  can  translate  and 
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rotate  in  all  directions  \inder  the  influence  of  the  external  forces  (i.e., 
gravity,  aerodynamic,  impact),  as  well  as  the  constraining  internal  element 
forces.  The  movement  of  the  masses  results  in  changes  in  the  relative  dis- 
tortion of  the  structural  elements  and,  in  turn,  results  in  a new  set  of 
element  forces  acting  throughout  the  system. 

Computer  Program  KRASH  has  the  capability  to: 

1.  Define  the  response  of  six  degrees  of  freedom  (DOF)  at  each 
representative  location,  including  three  translations  and 
three  rotations 

2.  Determine  mass  accelerations,  velocities,  and  displacements 

and  internal  member  loads  and  deformations  at  each  time  interval 

3.  Provide  for  general  nonlinear  stiffness  properties  in  the  plastic 
regime,  including  different  types  of  load-limiting  devices,  and 
determine  the  amovint  of  permanent  deformation 

U.  Determine  how  and  when  rupture  of  an  element  takes  place  and 
redistribute  the  loeuiing  over  the  structural  elements 
involved 

5.  Define  mass  penetration  into  an  occupiable  volume 

6.  Define  the  volume  change  due  to  structural  deformations  of  the 
occupiable  volume 

7.  Provide  for  ground  contact  by  external  structure  including 
sliding  friction  and  a nonrigid  ground  surface 

8.  Include  internal  structural  damping 

9.  Include  a measure  of  injury  potential  to  the  occupants;  for 
instance,  the  probability  of  spinal  injury  indicated  by  the 
Dynamic  Response  Index  (DRI) 

10.  Determine  the  distribution  of  kinetic  and  potential  energy  by 
mass  item,  the  distribution  of  strain  and  damping  energy  by 
element,  and  the  crushing  and  sliding  friction  energy  associated 
with  each  external  spring 

11.  Determine  the  vehicle  response  to  an  initial  condition  that 
includes  linear  and  angular  velocity  about  three  axes  and  any 
arbitrary  vehicle  attitude 
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12.  Provide  a measure  of  the  airplane  eg  velocity  by  means  of 
translational  momentum  relationships. 

13.  Analyze  an  impact  into  a horizontal  ground  and/or  an  inclined  slope. 

1^.  Provide  a measure  of  the  internal  stress  state  of  internal 
beam  elements 

15.  Analyze  a mathematical  model  containing  up  to  8o  masses  and 
100  internal  beam  elements 


1.3  THEORY 


1.3.1  Mathematical  Model  of  Airplane 

The  airplane  is  modeled  as  a series  of  interconnected  lumped  masses. 
Each  mass  is  allowed  six  degrees  of  freedom,  three  translations  and  three 
rotations.  The  masses  are  connected  internally  by  nonlinear  beam  elements. 
Each  mass  is  allowed  up  to  three  external  nonlinear  springs,  which  radiate 
outward  from  the  mass  and  contact  the  ground,  providing  external  crash 
forces.  A typical  airplane  model  is  shown  in  Figure  1-1. 


The  three  lumped  masses,  m^,  m^  and  m^^,  are  interconnected  by  beams 
ij  and  Jk.  Three  external  springs  are  shown  radia\.ing  outward  from  nij^.  The 
end  points  of  these  springs  are  denoted  by  C , , C and  C^,_.  The  center-of- 
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gravity  of  the  entire  airplane  is  denoted  by  G. 
no  lumped  mass  exactly  at  G.) 


"yk  “■  "zk" 

(In  general,  there  will  be 


1.3.2  Coordinate  Systems 

The  following  coordinate  systems  are  established  to  facilitate  the  deri- 
vation of  equations  for  the  mathematical  model. 

1.3. 2.1  Ground  Coordinate  System  (Oxyz) 

This  is  a right-handed  coordinate  system  fixed  in  the  ground  with  origin 
at  point  0 in  Figure  1-1.  The  x axis  is  positive  forward,  the  y axis  posi- 
tive to  the  right,  and  the  z axis  positive  downward.  The  xy  plane  (z  = O) 
corresponds  to  the  ground  surface . The  ground  coordinate  system  is  con- 
sidered an  inertial  coordinate  system  for  writing  the  dynamic  equations  of 
motion. 
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Figure  ] 


1.3. 2.2  Slope  Coordinate  System  {0x„  y z ) 

1 all  DL 

This  is  a right  handed  coordinate  system  fixed  in  the  ground  with  origin 
at  point  0 as  shown  in  Figure  1-1.  The  axis  is  positive  forward  up  the 
slope,  the  yg^  axis  is  positive  to  the  right,  and  the  Zgj^  axis  is  positive 
perpendicular  to  the  slope,  going  into  the  slope.  This  coordinate  system  is 
the  same  as  the  ground  coordinate  system,  rotated  through  an  angle  beta  (3)  posi- 
tive clockwise  about  the  Oy  gro\ind  axis.  The  0x„  y plane  represents  a 

bli  oL 

plane  inclined  at  an  angle  beta  with  respect  to  the  horizontal  ground 
plane.  Beta  is  a constant  input  angle  that  can  lie  in  the  range  from 
zero  to  ninety  degrees. 

1.3. 2. 3 Airplane  Coordinate  System  (Hx"y"z") 

This  is  a left-handed  coordinate  system  fixed  in  the  airplane  with  ori- 
gin at  point  H in  Figure  1-1.  The  x"  axis  is  positive  aft,  y"  is  positive 
left,  and  z"  is  positive  upward.  The  origin,  point  H,  corresponds  to  F.S.O, 
B.L.O,  W.L.O.  This  coordinate  system  is  used  only  to  input  the  coordinates 
of  the  lumped  mass  points,  since  the  coordinates  of  the  points  are  usually 
available  in  terms  of  F.S.,  B.L.  and  W.L. 

1. 3.2.1  Center-of-Gravity  Coordinate  System  (Gx'y'z’) 

This  is  a right-handed  system  fixed  in  the  airplane  with  the  origin  at 
the  vehicle  eg  (Point  G).  The  x'  axis  is  positive  forward,  y'  positive  right, 
and  z'  positive  downward.  These  axes  are  parallel  to  the  Hx"y"z’'  axes. 

1.3. 2.5.  Body  Coordinate  System  ) 

Each  lumped  mass  has  its  own  right-handed  coordinate  system  fixed  in 
the  mass.  The  initial  orientation  of  each  of  these  coordinate  systems  is 
arbitrary,  and  it  is  specified  by  means  of  three  input  Euler  angles  for  each 
mass,  relating  its  initial  orientation  to  the  Gx'y'z'  center-of-gravity  coor- 
dinate system.  Normally  the  body  coordinate  system  is  taken  as  initially 
parallel  to  the  center-of-gravity  coordinate  system  (and  hence,  the  three 
input  Euler  angles  are  set  =0),  since  the  inertia  data  is  generally  available 
about  these  eoces. 
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The  body  coordinate  system  is  the  system  used  in  writing  Euler’s  equations 
of  motion  for  each  lumped  mass  m^. 

1.3. 2. 6 Beam  Coordinate  System  ^ 


This  is  a right-handed  coordinate  system  with  the  xb  axis  along  a 

^ J 

straight  line  from  m.  to  m . As  the  masses  move,  the  coordinate  system 

^ J 

changes  orientation  so  that  xb  is  always  pointing  from  m.  to  m . If  the 

1 J 1 J 

beam  connects  node  points  which  are  offset  from  the  mass  points,  then  xb 

^ J 

always  points  from  the  node  point  rigidly  attached  to  mass  m^  to  the  node 
point  rigidly  attached  to  mass  m^ . 

The  direction  of  yb  and  zb  (they  are  mutually  perpendicular)  is 

1 J J.  J 

arbitrary,  and  is  defined  internally  within  the  program.  Each  beam  ij  has 
a beam  coordinate  system,  used  to  compute  the  beam  forces  and  moments. 

1.3.3  Relations  Between  Coordinate  Systems 

Any  quantity  which  can  be  expressed  as  a vector  in  one  coordinate  system 
(forces,  moments,  displacements,  velocities,  accelerations,  but  not  large 
rotations)  can  be,  likewise,  specified  in  another  coordinate  system  by  use 
of  an  Eulerian  transformation  matrix  relating  the  two  coordinate  systems. 
Thus,  for  example,  if  we  have  a force  vector  at  m^  expressed  in  its  body 
axes  components  as  (X^,  Z^^),  this  can  be  expressed  in  ground  axes  as 

simply 


V Oi  / Vi/  (1_1) 

where  1 X . , Y . , Z , I are  the  ground  axes  components  of  the  force  vector  and 
1 oi  oi  oi  » - 

Ia^I  is  a 3 X 3 Eulerian  transformation  matrix.  The  form  of  lAj  depends 
on  the  sequence  in  which  rotations  from  one  axis  system  to  the  other  are 
performed.  In  this  program,  the  following  sequence  is  employed: 

1.  Starting  with  the  Oxyz  coordinate  system,  perform  a clockwise 

(right-hand  rule)  rotation  about  the  Oz  axis,  through  an  angle  4^ 
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2.  Next,  perform  a clockwise  rotation  about  the  new  Oy  axis,  through 
em  angle 

3.  Finally,  perform  a clockwise  rotation  about  the  final  Ox  axis, 
through  an  angle 

The  resulting  transformation  matrix  is  then  inverted  to  obtain  a trans- 
formation from  body  axes  to  groiind  axes.  The  resulting  trems format ion  matrix 
is  given  by 


COS©^  C0S^^|-C0S^^  +SIN^5^  SINGj^  COS^jj  SIN^5^SIN^j+COS^5^SIN0^COS^j 

COS©^  SiN^^I  COS0^  +SINj^^  SIN©^  SIN^J-SIN?5^C0S«|»^+C0S{ijSIN©^SIWi/f^ 


-SIN©, 


SIN^J^  COS©^ 


COS^fj^  COS©^ 


(1-2) 


Similar  transformation  matrices  are  formed  to  relate  various  other  coor- 
dinate systems.  These  are  summarized  in  Table  1-1. 


TABLE  1-1.  TRANSFORMATION  MATRICES 


Matrix 

Transforms  From 

To 

Using  Angles 

Angles 
Constant  ( 
Varying 

IBB 

ith  body  axes 

ground  axes 

^1» 

varying 

[*u] 

beam  ij  axes 

ground  axes 

•^ij*  ®ir  '^'ij 

varying 

[*:] 

ith  body  axes 

eg  axes 

constant 

[*•] 
r T 

eg  axes 

ground  axes 

«!>',  6', 

constant 

N 


ground  axes 


slope  axes 


0,  p , 0 


constant 


r 


All  the  above  matrices  utilize  equation  (l-l)  with  the  appropriate 
angles  from  Table  1-1.  and  |a'J  are  used  only  in  initial  conditions 

determination,  so  that  their  Euler  eingles  are  Indicated  as  constant.  These 
"constants"  are  really  the  time  zero  values  of  the  angles  which  actually  vary 
but  are  not  used  in  the  time  history  computations.  <f>' , 0',  i|)'  are  input 
constants  defining  the  initial  attitude  of  the  overall  vehicle. 

We  also  will  require  the  time  derivative  of  the  [\]  matrix,  denoted 
It  can  be  shown  that  this  is  given  by  postmult iplying  ^ by  a 


[M- 

matrix 


so  that 


N 

[*i]  - [*i][“i] 


(1-3) 


where 


[“0 


is  the  following  skew-symmetric  matrix: 


N 


-^i  ^i 


-p^ 


-q.  p. 

1 X 


(l-U) 


P^t  are  the  body  axis  components  of  the  angular  velocity  of  mass  m^. 

1.3.**  Sign  Conventions 

The  basic  sign  convention  used  for  all  displacements,  rotations,  veloc- 
ities, accelerations,  forces  and  moments  is  that  all  quantities  are  positive 
in  the  positive  direction  of  the  axes  shown  in  Figure  1-1.  Rotation  and 
moments  utilize  a right-hand  rxile  to  define  the  positive  direction.  The  few 
exceptions  to  this  rule  are  indicated  in  the  theoretical  development . 
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1.3.5  Forces  Acting  on  Each  Mass 


The  fol.lovinp:  forces  and  moments  act  on  each  lumped  mass: 

• Gravity  forces 

• Aerodynamic  forces 

• Internal  forces  and  moments  - Strain  and  Damping 

• External  forces  and  moments 

All  forces  and  moments  are  positive  when  they  act  in  the  positive 
direction  of  the  axes. 


1.3. 5.1  Gravity  Forces 

The  gravity  force  for  the  ith  lumped  mass  is  Just  the  weight  acting 
along  the  ground  fixed  Oz  axis.  Transforming  this  force  into  body  axes 
gives  the  gravity  forces  as 


/ X.  ' 

I 

0 \ 

1 ' 

1 1 

) 

I Y 

1 0 ( 

1 I 

j 1 

\ ^ ] 

where  indicates  the  transpose  of  Note  that  the  transpose  rather 

than  the  inverse  is  used,  since  for  rotation  transformation  matrices  the  two 
are  equal. 

1.3. 5.2  Aerodynamic  Forces 

The  only  aerodynamic  forces  considered  are  a constant  lift  force,  in 
ground  axes,  positive  upward,  acting  on  each  lumped  mass.  This  lift  force 
is  expressed  as  a constant  fraction  of  the  total  vehicle  weight,  times  the 
vehicle  weight : 


LIFT^  = 1 

i TOT 


(l-o) 


where 


'^TOT  “ '^i 


(1-7) 


and  N is  the  total  number  of  lumped  masses.  Transforming  this  lift  into  ith 
body  axes  gives  the  desired  aerodynamic  forces  as 


(1-8) 


- LIFT. 


where  the  minus  sign  is  necessary  because  the  lift  is  positive  upward  (in 
the  negative  direction  of  the  Oz  ground  axis). 

1 . ■( . . J Internal  Forces  and  Moments  - Strain  and  Damping 

These  forces  and  moments  result  from  the  deformations  of  the  inter- 
connecting beams  as  the  lumped  masses  move.  First,  for  beam  i j , the  deflec- 
tions, rotations  and  velocities  of  beam  end  points  j and  i are  computed  in  beam 
ij  axes.  These  are  then  used  to  calculate  the  three  forces  and  three  moments 
it  I'oint  ,i , utilizing  equations  based  upon  the  12  x 12  stiffness  matrix  relat- 
ii.,  t'.ie  loads  t(  i;e  deformations.  Nonlinear  effects  are  also  included  in 
co::;nutinf  the  loads.  Tie  loads  at  point  i in  beam  axes  are  computed  from 
static  equi  1 ibriuiii  equations  (the  beams  are  assumed  massless).  Then  the  loads 
:.t  i.  and  J are  transfornu'.l  into  ith  and  jth  mass  body  axes,  respectively. 


1-10 


Finally,  for  each  mass  the  contributions  from  the  various  beams  that  attach 
to  that  mass  are  summed  to  obtain  the  total  internal  forces  for  the  mass. 

All  the  internal  'aln  force  calculations  are  performed  on  an  incremental 
basis,  i.e.,  incremental  deflections  are  used  to  compute  incremental  forces, 
since  a tangent  modulus  is  utilized  for  computing  nonlinear  effects.  The 
damping  forces  are  total,  not  incremental,  based  on  current  beam  end  point 
velocities . 


1.3. 5.3.1  Beam  Deformations  and  Velocities 


The  first  step  in  determining  the  internal  beam  strain  and  damping  loads 
is  to  calculate  the  deflections  and  rotations  of  the  end  points  of  the  beam, 
as  well  as  the  velocities  of  these  points.  Figure  1-2  illustrates  a typical 
beam  whose  position  in  space  is  defined  by  the  ground  coordinates  of  the  end 
points  of  the  beam.  This  figure  illustrates  the  most  general  case  in  which 
the  beam  connects  two  node  points  which  are,  in  turn,  rigidly  connected  to 
mass  points  i euid  j.  Beams  can  also  connect  directly  to  the  mass  points. 

The  purpose  of  the  node  points  is  to  allow  a finite  size  for  the  mass  points 
so  that  beam  end  connections  can  be  offset  from  the  lumped  mass  center  of 
gravity. 


Masses  i and  j are  located  at  (x.,y.,z.)  and  (x  ,y  ,z  ) respectively, 

111  J J J 


where  these  are  the  ground  axes  coordinates  of  m.  and  m 

^ J 


These  are  time- 

varying  quantities  obtained  from  numerical  integration  of  the  equations  of 
motion  for  m.  and  m 


J 


Node  points  iM  eind  jN  are  located  at  (xnp^, 


and  (xnp  , ynp,„,  znp,^),  respectively.  These  coordinates  are  calculated  as 
J 


follows 


. 

1 xnp 

X 

rx 

jynp 

= 

y 

- [AJ 

ry 

I 2np 

iM 

z 

i 

rz 

(1-9) 


[A^]  is  given  by  Equation  (1-2)  and  is  a time  varying  quantity  dependent 

upon  the  Euler  angles  0^ , 0^,  emd  that  define  the  orientation  of  mass  m^ 

will)  respect  to  ground.  The  vector  (rx,ry,rz)  is  a constant  defining  the 

iM 
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o:‘:'se:  o:'  r.cde  ycir.t  !;•;  with  respect  to  mass  m^;  these  vector  components  are 
in  mass-:'ixed  axes. 

Taking  the  time  derivative  of  both  sides  of  Equation  (1-9)  yields 


r ^ N M r 


(1-9) 


where  Equation  (l-3)  has  been  used  for  • Multiplying  through  by  the 

transpose  of  convert  from  ground  axes  back  to  mass  axes  yields 


(1-10) 


'•■"P  t iM 


The  (u,v,w)^  and  (unp,vnp,wnp)^j^  are  the  body  eixes  components  of  the 
velocity  vectors  for  mass  i and  node  point  iM,  respectively.  The  former  is 
obtained  from  integration  of  the  equations  of  motion. 

Since  the  strain  forces  eire  calculated  on  an  incremental  basis,  what  is 
required  are  the  incremental  deflections  of  the  end  points  of  the  beam.  These 
are  calculated  directly  as  the  differences  between  the  current  coordinates 
£Uid  the  values  at  the  previous  time  point. 


(1-lla) 


1 dx 

1 

X 

1 dy 

= 

y 

- 

1 dz 

i 

z 

i 

current 


previous 


dxnp 

xnp 

1 

xnp 

dynp  = 

ynp 

- 

ynp 

iM 

znp 

1 iM 

znp 

d-llb] 


current 


\ I iM 
previous 
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The  incremental  rotations  are  also  treated  as  vectors.  Although  large 
rotations  are  not  vector  quantities,  the  incremental  rotations  from  one  time 
step  to  the  next  are  sufficiently  small  to  be  treated  as  vectors.  The  angular 
velocity  components  of  m^  are  defined  as  {p,q,r)^  in  mass  fixed  axes.  These 
are  obtained  from  integration  of  the  angular  accelerations  {p,4»r)^*  The 
incremental  rotations  are  calculated  as  the  differences  (ctirrent-previous) 
in  the  integrals  of  the  angular  velocities. 


dpin 

inp 

dqin 

inq 

drin 

i 

inr 

i 

current 


inp 
inq 
inr 

previous 
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where 


inp 

p 1 

inq 

= f 

q I 

inr 

1 J 

' 1 o 

r 1 

dt 


(1-13) 


The  incremental  rotations  in  Equation  (l-12)  are  in  mass-fixed  axes  and 
follow  the  right  hand  rule. 

In  the  actual  computer  program  KRASH,  the  internal  beam  strain  and  damping 
forces  are  calculated  in  a DO  loop  ranging  over  eJ.1  the  beams  in  the  model 
being  analyzed.  For  each  beam,  the  end  points  eire  defined  by  mass  numbers  i 
and  J amd  node  point  numbers  M and  N.  If  M or  N is  zero,  the  beam  connects 
directly  to  the  mass  point  at  that  end.  The  appropriate  equations  are  then 
used  to  calculate  the  end  point  positions,  velocities  and  incremental,  chauiges 
in  position.  For  the  remaining  equations,  the  following  terminology  is  employed: 


dxmod 

j 

dx 

dxnp 

dymod 

= 

dy 

1 or  ( 

dynp 

dzmod 

1 

i 

dz 

k 1 

dznp 

iM 


(1-lUb) 
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umod 
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u J 

unp 

vmod 

V 1 or  ( 

vnp 

wmod 

li  1 

’'li  1 

vnp 

iM 


(1-ll+c) 


In  each  case,  the  first  terms  are  used  if  there  is  no  node  point  (M  = 0 
for  the  current  beam),  and  the  second  terms  are  used  if  there  is  a node  point 
(M  ^ 0).  For  the  angular  velocities  (p,q,r)j^  and  the  incremented,  rotations 
(dpin,  dqin,  drin)^,  the  quantities  are  the  same  whether  or  not  there  is  a node 
point,  since  the  node  point  iM  is  rigidly  connected  to  mass  i.  Beam  end  j or 
jN  is  treated  in  the  same  fashion. 

Once  the  quantities  in  Equations  (l-lU)  have  been  calculated,  the  beam 
orientation  with  respect  to  ground  is  known  so  that  the  current  [a  1 matrix 
may  be  determined.  Figure  1-3  illustrates  a typical  beam  i j , showing  its 
spatial  orientation  with  respect  to  the  ground  fixed  coordinate  system  and 
identifying  the  beam  Euler  euigles  0 and  From  Figure  1-3,  it  can  be 

seen  that 


xiJ 

i xmod 

xmod 

yiJ 

= 1 ymod 

ymod 

ziJ 

1 1 zmod 

1 J 

zmod 

(1-15) 


ij 


tan“^  (yij/xij) 


{l-l6a) 


®ij  = - tan"^  (zij/ 


5 2 

xij  + yij 


(l-l6b) 


The  remaining  Euler  angle  ^ , defining  the  "roll"  orientation  of  the 

^ J 

betun  axes,  is  initially  set  to  zero  and  thereafter  calculated  as 
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(l-l6c) 


current 


d + dd 
'"ij 

previous 


where  dd.,  is  the  average  incremental  rotation  of  masses  i and  j resolved  into 
the  current  beam  x axis.  This  term  is  derived  later,  and  is  given  by  Equation 
(l-3d).  The  current  beam  length  is  calculated  as 


£ = >/ 

^.2  ..2  2 
xij  + yij  + zij 

(1-17) 

In  order  to  transform  vector  quantities  directly  from  mass 

beam  axes,  the  following  transformation  matrices  are  calculated 

fixed  axes  to 

[aijtai] 

(l-l8a) 

[aijtaj] 

(l-l8b) 

In  Equations  (l-l8),  the  second  matrix  transforms  from  mass  fixed  axes  (either 
i or  j ) to  ground  axes,  and  the  first  matrix  transforms  from  ground  axes  to 
beam  ij  axes. 

The  incremental  beam  deflections  at  i «ind  J,  in  beam  iJ  axes,  are  given  by 


dx 

dy 


I 


dxmod 

dymod 

dzmod 


The  beam  end  translational  velocities  are  simply 


vx 

umod 

vy 

= [aijtai] 

vmod 

vz 

i 

wmod 

(1-19) 


(l-20a) 
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The  incrementatl  rotations  of  the  beam  ends  are  calculated  as 
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The  beam  end  rotational 

velocities  are 
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(l-21a) 


(l-21b) 


(l-22a) 


(l-22b) 


1.3.5«3.2  Linear  Strain  Forces  and  Moments 

The  incremental  beam  strain  forces  and  moments  are  first  calculated  on 
a linear  basis  and  then  modified  by  a stiffness  reduction  factor  KR,  which 
may  vary  at  each  time  step.  The  calculation  of  the  KR  factors  is  discussed 
in  subsection  1.3. 5. 3. 4;  this  section  presents  the  method  of  calculating  the 
linear  (KR=1)  strain  forces  and  moments. 


1-1 8 


The  corresponding  incremental  force/deflection  equation  is 


';1  ■ l~;! 


(1-2U) 


where 


{dVi)  S 


1 dx' 
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dx\ 

1 dy 

idy  1 

Idz 
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dz  1 

|d4.j 

1 de 

1 j 

del 

i 

(1-25) 


The  terms  in  Equations  (1-25)  are  given  by  Equations  (1-19)  and  (l-2l). 
|dFi|  and  |dFj}  are  6 element  vectors  for  the  forces  and  moments  at  beam  ends 
i and  J;  these  forces  are  ordered  in  the  same  manner  as  the  deflections,  i.e. 


(^4  ' 


(1-26) 


All  quantities  in  these  equations  are  in  beam  ij  axes. 

If  Equation  (1-2U)  is  multiplied  out  using  (1-23),  the  axial  and  torsional 
(x  and  «<> ) relationships  are  simply 


dM<|^  = (d4)j-d4)^)  , dM<l5^  = -dM4>j 


(l-27a) 

(l-2Tb) 
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(l-28a) 

(l-28b) 


^ AE/e 

GJ/l 
~ > 

The  other  forces  and  moments  consist  of  two  sets  of  coupled  bending  equa- 
tions. Deflections  in  the  z direction  and  0 rotations  are  coupled  in  forming 
dFz  and  dM6,  and  y deflections  and  4^  rotations  are  coupled  in  the  equations 
for  dpy  and  dM4^.  The  following  relationships  result  from  expanding 
Equation  (l-2U): 


3kQ 

dMBj  = k20(dZj  - dzj^)  +—  (dSj  - d0j^)  + — (dB^  + de^ ) (l-29a) 

dFz^  = “dFZj 

dMBj^  = -dMGj  + ^dFZj  (l-29b) 

where 


k,  = 12EIy/g^ 
z 

= 6EIy/^2 

kg  = UEIy/£ 

The  corresponding  equations  for  the  coupled  y,  4^  motions  are 


(l-29c) 


•ilVj  = ^"^i^ 

dMVj  = lV4(^‘^J  " ‘^i^ 

dFy^  = -dFyj 

dM4<j^  = -dM4>j  - idFyj 


(l-30a) 


(l-30b) 
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where 


ky  = 12EIz/«^ 
ky^  = -6eIz/«^ 

k^  = kEIz/2  (l-30c) 

In  the  above  equations,  E is  the  modulus  of  elasticity,  G is  the  modulus 
of  rigidity,  £ is  the  current  beeun  length,  ly  and  Iz  are  the  beam  cross-sectional 
area  moments  of  inertia  about  the  beam  y and  z axes,  respectively,  and  A is 
the  cross-sectional  area.  J is  the  torsional  stiffness  factor;  for  circular 
cross  sections  this  is  equal  to  the  polar  moment  of  inertia  (ly  + Iz),  while 
for  non-circular  sections  J is  less  than  the  polar  moment  of  inei^ia. 

In  program  KRASH,  the  following  incremental  deflection  sums  and  differ- 
ences are  defined 


(l-31a) 


(l-31b) 
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utilizing  these  definitions,  the  incremental  force  and  moment  equations 
become 


dFx,  = k d, 
J X 1 


dFVj  = k^  d^  . ky^d^ 
dFZj  - k^  d3  + dg 

= k^  dj^ 


dMe.  = k^e  S ^5  H 


While  the  corresponding  quantities  at  the  i end  of  the  beeun  are 


dFXj^  = -dFXj 

dFy^  = -dFyj 

dFz^  = 
dM4>^  = -dM4>j 

dMe^  = -dMGj  + idFZj 

dM4>^  = -dM4<j  - £dFy-j 


(l-32a) 


(l-32b) 


Figure  1-U  Illustrates  these  incremental  forces  euid  moments  at  both  ends 
of  beam  i j ; positive  values  are  shown  in  the  figure. 

In  Equation  (l-3l),  the  d^  term  is  used  to  update  the  torsional  orientation 
of  the  beam  axis  system.  Equation  (l-l6c)  indicates  that  the  beam  Euler  angle 
4>. , is  continuously  updated  by  <i<^.  ,•  This  is  taken  as  the  overage  incremental 
rotation  of  the  1 and  J ends  of  the  becun,  and  is  given  by 


d+ij  = d,/2 


Where  d^  is  given  by  Equation  (1-31)  as  (d4>j  + d<t>^) 
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(1-33) 


r 


[ Figure  1-4. 


Beam  Horizontal  Strain  Forces  emd  Moments 
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1.3.5«3.3  Pinned  End  Conditions 


Equations  (l-32)  developed  in  the  preceding  section  are  applicable  only 
for  the  condition  wherein  the  beam  is  built  in  or  fixed  at  both  ends.  If  the 
beam  is  pinned  at  both  ends,  then  both  the  force  and  moment  at  1 and  J are  zero 
in  the  plane  perpendicular  to  the  pin  axis.  If  the  beam  is  pinned  at  one  end 
euid  fixed  at  the  other  end,  then  the  equations  derived  in  the  preceding  section 
are  not  applicable. 

In  order  to  illustrate  the  derivation  of  the  force-deflection  relation- 
ship for  a single  pinned-end  condition,  the  coupled  z - 0 bending  equations 
will  be  considered.  Isolating  these  terms  from  Equation  (l-23),  we  have 


r \2H^  -611^ 


El 

y 


kt2 

( symmetric ) 


-12/i^  -6/11^ 
6tt  ^ 2li 
12111^  6/11^ 
h/i  _ 


(1-3U) 


If  there  is  a pinned  connection  at  J,  then  dM6j  is  zero,  and  the  above 
equation  can  be  rewritten  as 


= El 

^.^21^^22.  '^"'^2 


(1-35) 


Kii  dv^  + K^2  <^^2  " 


^21  *^^1  *^22  *^^2  “ ° 


Solving  the  second  equation  for  dVg  and  inserting  into  the  first  equation  yields 

[■Si  - "12  V""  "21]  '*''1 ' 


12  22  21 


‘1'  y 
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Rearranging,  ve  have 


where 


= El  [ic]  dv^ 


(1-37) 


dF,  t < dMe 


dvi  4 U 


(Kl  4 [lC,J  - [k^,]  [kJ-^  [kJ 

■ 12/P  - 6/P  - 12/p' 

* -6/£^  k/t  6/P  (3x3  matrix) 

.-12/P  6/P  12/P. 


(l-38a) 


(l-38b) 


(l-39a) 


N = \ 2/2  / (3  X 1 column  matrix) 

( 6/P  ) 

f^2lj  ~ 6/2^]  (1x3  row  matrix) 

[^^22j~  (1x1  scalar) 

Incorporating  (1-39)  into  (1-38),  (1-37)  becomes 


■ 1/2^  -1/P  - I/Pt  JdZj 

3EIy  -1/P  1/2  1/2^  <dej 

-1/P  1/P  1/2^  Idz, 


(l-39b) 


(l-39c) 


(l-39d) 


(I-UO) 


utilizing  the  stiffness  teras  defined  by  Equation  (l-29c)  this  reduces  to 


T ’'e 


Multiplying  this  olit  and  using  the  definitions  of  (l-3l),  we  obtain 


dFzj 


dM0j  = 0 

dFz^  = -dFZj 


dMG^ 


d0. 


( Pinned  at  J ) 


(1-1+la) 


Similar  derivations  yield  the  following  equations  for  the  situation  where  the 
beam  is  pinned  at  i and  fixed  at  J. 


(Pinned  at  i) 

(l-l*lb) 

dFz^  = 

dM0^  = 0 

Equations  (l~l*l)  are  also  applicable  for  the  coupled  y - ^ bending  if  subscripts 
are  changed  appropriately.  The  resulting  equations  are  shown  for  reference: 
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■ "f  ^2  * ■f*' 

“+1  - 2 ^2  * I ’‘t 

dFy^  = -dPVj 
dM4p^  = 0 

In  program  KRASH,  only  the  equations  for  the  forces  and  moments  at  J are  calcu- 
lated directly  from  (1-1*1)  and  (1-1*2).  The  forces  and  moments  at  i are  calcu- 
i 'w'  lated  using  the  static  beHance  Equations  (l-32b);  it  ceui  be  shown  that  the 

results  are  the  same  as  the  i end  equations  from  (l-l*l)  and  (1-1*2). 

This  concludes  the  discussion  of  the  linear  incremental  strain  forces  and 
< moments;  the  following  section  describes  the  treatment  of  the  nonlinear  behavior 

of  the  beam  elements. 

( 

1.3.5. 3.1*  Nonlinear  Strain  Forces  and  Moments 

Figure  1-5  illustrates  a typical  nonlinear  load-deflection  characteristic 

» 

curve,  in  this  case  the  axial  load  at  point  J,  Fx j » versus  the  axial  deflection 
of  J with  respect  to  i,  vx.  Program  KRASH  utilizes  a piecewise  linear  repre- 
sentation of  the  load-deflection  curve.  At  each  increment  in  time,  a linear 

I 

incremental  force  <iF'Xj  is  calculated  according  to  the  equations  of  the  previous 
sections.  This  incremental  force  is  then  multiplied  by  a stiffness  reduction 
factor,  KR  , to  obtain  the  nonlinear  incremental  force  dP'  . These  increments 

X , X 

are  then  summed  over  time  to  form  the  total  force  Fx^.  Figure  1-5  illustrates 
a typical  incremental  deflection  dvx,  the  corresponding  linear  incremented 
force  dFXj  and  the  final  nonlinear  incremental  force  dFxj  . 
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The  stiffness  reduction  factor  KR^  is  an  input  tabular  function  of  deflec- 
tion vx.  Figure  1-6  shows  the  KR^  vs  vx  curve  that  corresponds  to  the  load- 
deflection  curve  in  Figure  1-5.  Clearly  KR^  is  just  the  first  derivative  of 
the  load-deflection  curve  with  respect  to  deflection.  Any  KR  curve  can  be 
input  by  the  program  user  to  model  any  general  nonlinear  behavior.  Certain 
standard  curves  are  available  to  represent  the  more  common  types  of  nonlineari- 
ties. These  are  discussed  in  detail  in  the  User's  Guide  (Volume  II). 

Separate  KR  tables  are  input  for  each  of  the  6 directions  for  which  it  is 
desired  to  model  nonlinear  behavior,  for  beam  i j . In  practice,  perhaps  only 
20X  or  less  of  all  the  beams  will  be  modeled  nonlinear ly.  (If  no  KR  table  is 
input,  KR  = 1 and  the  beam  is  linear  for  the  entire  run.)  In  general  then,  we 
have 


where  the  diagonal  [KR]  matrix  is  composed  of  six  stiffness  reduction  factors, 
each  of  which  is  an  input  tabular  function  of  the  corresponding  total  beam 
deflection.  The  KR's  are  functions  of  the  time  sums  of  the  first  six  incre- 
mental deflections  in  Equation  (l-3l);  these  are  all  the  (j  - i)  relative 
deflections ; 
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REDUCTION 
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KR« 


/ i. 


d«x  — 


dFx  - KRx,  dFxj 


AXIAL  DEFLECTION  *x 


Figure  1-5.  Nonlinear  Load  Deflection  Curve 


,KR„  • 1.0  CORRESPONDS  TO  LINEAR  REGION 


y KRx  - 0 CORRESPONDS  TO  PLASTIC  REGION 


= f(vx) 

KRy.  = f(vy) 

KR  = f(vz) 
z 

KR^  = f(v4>) 
krq  = f(ve) 
KR^  = f(v4>) 


(1-U5) 


The  toteLL  beam  forces  are  calculated  as  the  time  sums  of  the  incremented,  forces 
from  Equation  (I-U3) : 

t 

(l-U6a) 

I *J  ; 
o 


' i h'} 


{V}  = I {“'1') 


(1-U6b) 


These  total  strain  forces  are  in  beeira  axes  and  act  at  the  ends  of  the  beeun 
as  shown  in  Figure  1-^*  for  the  incremented,  strain  forces. 

For  the  uncoupled  degrees  of  freedom  (x  and  4) ) , the  force  depends  only  on 
one  relative  deflection  and  the  corresponding  KR  curve.  This  situation  was 
illustrated  previously  for  the  x (axial)  direction.  For  the  coupled  degrees  of 
freedom.  Equation  (l-32a)  implies  that  a given  force  depends  on  2 or  3 deflection 
terms  and  the  one  KR  cujrve  corresponding  to  the  force.  For  example,  if  we  com- 
bine Equations  (l-32a)  and  (1-43)  for  the  coupled  z,  6 degrees  of  freedom, 
we  have 


apzj'  . [k^  dj  . k^e  dg] 

k 

dH6j'  = KR^gCdg)  [k^g  dg  ♦ aj 


(1-47) 


where  the  parentheses  after  the  KR  terms  indicate  their  functional  relationships , 
not  a multiplication. 


1 
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If  the  calculations  are  actually  performed  in  the  above  manner,  the  inde- 
pendent behavior  of  dz  and  d6  (and  therefore  d^  and  d^)  will  result  in  differ- 
ent values  for  and  KR0  at  any  point  in  time.  Since  the  proportionality  of 
dFz  and  dM6  that  always  exists  in  the  linear  system  is  lost,  this  computational 
scheme  has  a strong  tendency  to  ;d.eld  negative  strain  energy,  which  a real 
physical  system,  even  nonlinear,  can  never  do.  Furthermore,  one  would  intui- 
tively expect  that  any  failure  mechanism  in  the  actual  structure  which  causes 
yielding  or  nonlinear  behavior  in  the  z direction-  would  also  cause  similar  non- 
linear behavior  in  the  coupled  6 direction.  This  reasoning  could  be  extended 
to  include  all  six  independent  deflections,  so  that  when  any  one  direction  dis- 
plays nonlinear  behavior,  this  same  nonlinearity  holds  for  all  6 directions. 
However,  it  is  felt  that  the  argument  for  similar  behavior  is  most  persuasive 
for  the  coupled  degrees  of  freedom,  and  that  the  axial  and  torsional  modes 
should  be  treated  independently. 

In  light  of  the  foregoing,  the  program  is  coded  in  the  following  manner. 
For  the  coupled  degrees  of  freedom,  the  direction  which  goes  nonlinear  first 
determines  the  nonlinear  behavior  for  both  coupled  directions.  In  the  above 
example,  this  means  that  if  KR^  departs  from  1 (linear)  before  KRg  does,  then 
is  used  for  both  drz  and  dM0'  for  the  remainder  of  the  run.  If  KR0  goes 
nonlinear  first,  then  it  is  used  in  calculating  both  dFz’  and  dM0'  for  the  rest 
of  the  run.  This  analytical  model  preserves  the  proportionality  that  exists 
between  dFz  and  dM0  in  the  linear  system,  since  the  linear  incremental  forces 
are  always  multiplied  by  the  same  KR  factor.  The  coupled  degrees  of  freedom 
y and  4>  are  treated  in  the  same  manner.  In  each  case,  both  deflection  and 
rotation  KR  curves  are  input,  even  though  only  one  will  be  used,  because  it  is 
not  known  in  advance  which  direction  will  predominate.  Axial  and  torsionsLL 
deflections  are  uncoupled  and  each  has  its  own  KR  curve,  applicable  for  the 
entire  analysis. 

In  addition  to  the  nonlinear  load  capability  Just  described,  the  program 
also  allows  for  unloading  and  subsequent  reloading  to  proceed  along  an  elastic 
line.  This  is  illustrated  in  Figure  1-7.  Th_  ordinate  Fxj*  is  the  total 
(over  time)  axial  load  for  beam  ij , from  Equation  (l-U6).  Similarly,  vx  is 
the  corresponding  total  (over  time)  axial  deflection  from  Equation  (l-UU). 
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Loading  proceeds  along  the  solid  line  0-F-A.  Unloading  then  proceeds 
along  dashed  line  A-B,  parallel  to  the  original  linear  loading  line  0-F, 
Subsequent  reloading  then  proceeds  along  B-A,  until  point  A is  reached.  At 
this  point,  further  loading  proceeds  along  A-C. 

Note  that  if  unloading  occurs  along  the  original  elastic  portion  of  the 
curve  0-F,  unloading  will  be  along  F-0.  The  intent  of  this  loading- 
unloading model  is  to  represent  crushing  structure  with  elastic  rebound. 

The  important  feature  represented  here  is  the  energy  absorbed  by  the  struc- 
ture in  a loading-unloading-reloading  cycle  such  as  0-F-A-B-A-C.  It  is  felt 
that  under  the  large  deformation  crash  conditions  that  this  prograin  is  intended 
to  investigate,  this  energy  absorption  is  far  more  significant  than  the  small 
structural  damping  generally  present. 

In  program  KRASH,  each  of  the  six  relative  (j  - i)  total  deflections 
(vyt  ...  viV)  is  monitored  for  unloading  and  reloading.  While  a given  deflection 
is  in  an  unloading  or  reloading  phase,  the  corresponding  KR  is  set  to  1.  This 
corresponds  to  the  original  linear  slope  as  with  A-B  in  Figure  1-7.  Subsequent 
to  this , in  order  to  preserve  the  linear  proportionality  for  the  coupled  de- 
grees of  freedom,  KRq  is  set  equal  to  KR^  and  KR^  is  set  equal  to  KR^,.  This 
latter  constraint  holds  true  even  if,  for  example,  0 is  unloading  and  z is 

loading  with  KR  j^l.  In  this  case,  KR  ?^1  will  dominate,  and  KRo  will  be  set 
z z ^ 

to  this  value  rather  than  1 even  though  0 is  unloading. 
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While  the  calculation  procedure  Just  described  does  not  guarantee  the 
absence  of  negative  strain  energy  for  nonlinear  coupled  degrees  of  freedom,  it 
has  been  found  by  experience  to  yield  acceptable  results  (positive  strain  energy) 
except  for  unusual  conditions  of  coupled  nonlinear  unloading  and  reloading.  It 
would  appear  that  the  stiffness  reduction  factor  approach,  which  is  the  heart 
of  program  KRASH's  philosophy  of  user  specified  element  load-deflection  be- 
havior, is  not  sufficiently  rigorous  to  eliminate  the  possibility  of  achieving 
negative  strain  energy  with  coupled  nonlinear  degrees  of  freedom.  A problem 
formulation  based  on  stress-strain  relationships  and  the  compatibility  of  ele- 
mental strains  and  beam  deflections  would  be  required.  Such  a formulation 
would  involve  considerably  more  computational  detail  than  is  used  in  program 
KRASH.  and  would  substantially  slow  down  the  program  execution.  In  light  of 


the  high  degree  of  success  with  practical  applications  of  program  KRA5H,  it  is 
felt  that  the  current  fcrnulation  is  quite  satisfactory  from  n practical  analy- 
sis standpoint,  if  not  from  a rigorous  analytical  viewpoint. 

Finally,  it  should  be  noted  that  only  the  coupled  nonlinear  degrees  of 
freedom  pose  potential  problems.  All  linear  motion  is  rigorously  correct  and 
never  yields  negative  .strain  energy,  and  the  uncoupled  degrees  of  freedom 
(x  and  4>  ) are  well-behaved  even  with  nonlinearities,  f'lnce  the  strain  energy 
of  each  beam  is  part  of  the  output,  the  user  will  be  aware  if  any  beams  with 
negative  strain  enerpy  exist  in  the  model. 

Program  KRAi'H  also  provides  for  complete  element  failure,  termed  rupture, 
if  any  of  the  six  beam  leliltive  deflections  {v{  or  six  total  beam  forces  jl'j'} 
exceed  specified  input  values.  Tf  any  deflection  or  force  exceeds  nllowable 

values,  that  i.l  beam  element  r\ipt\ires  and  all  its  forces  and  moments  are  set 
to  zero  for  the  remainder  of  the  run. 

l.l.S.t.i'  Unsymmetr leal  Beam  Klements 

The  befun  elements  described  previo\isly  bel\avc  symmetrically , meaning  tiiat 
the  linear  and  nonlinear  load-deflection  characteristics  are  the  sfune  for  botti 
j>oslt  Ive  and  negative  lieflecttons . Positive  deflections  yield  jioBitlve  loads, 
and  equal  negative  iJeflectlons  result  in  equal  negative  loads.  For  ceri.aln 
modeling  sitiiations,  it  may  be  desirable  to  have  available  be»im  elements  wltii 
different  charnctri'istlcn  for  ptu’.it.lve  and  negative  deflections.  While  pro- 
viding tills  capability  foi-  the  co\ipled  bending  degrees  of  freedom  would  be 
unduly  complicated  for  the  potential  benefits,  the  uncoupled  axial  degree  of 
fi'eedom  can  be  treated  in  a stralglit  forward  manner.  Tlie  axial  degree  of  free- 
liom  Is  also  the  one  wherein  \insymmel  rical  behavior  is  most,  conunnn. 

Situations  calling  for  unsymmetrlcal  axial  load-deflection  behavior  in- 
clude "string”  type  elements  such  as  seat  belts  or  wire  braces  and  diagonal 
elements  representative  of  shear  panel  behavior.  These  applications  call  for 
a tension-only  element,  with  zero  compression  capability.  Contact  between 
surfaces  not  originally  touching  can  be  modeled  with  compression-only  elements. 
An  example  of  thin  wo\»ld  be  an  element  connecting  a landing  gear  wheel  to  the 


( 


f; 

r 
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fuselage  or  vlng  structure.  Until  large  deformations  occur,  these  surfaces 
are  not  in  contact.  However,  after  the  failure  of  a landing  gear,  a wheel  may 
subsequently  contact  the  fuselage  or  wing.  A compression-only  element  with  a 
deadband  could  be  used  to  model  this  type  of  contact. 


Figure  1-8  illustrates  the  general  form  of  unsymmetrical  axial  load- 
deflection  characteristic  that  is  allowed  in  prograun  KRASH.  A tension-only 
element  with  backlash  is  illustrated.  The  load  is  zero  for  axial  extensions 
less  than  DBjj » an  input  deadbeuid  constant  for  beam  ij.  Beyond  this  point, 
the  normal  linear/nonlinear  beam  characteristics  are  employed,  with  the  follow- 
ing exception.  The  total  (over  time)  load  is  never  allowed  to  be  negative  (for 
a tension-only  element).  Therefore,  during  unloading  the  path  A-B-C  is  followed, 

holding  Fx  ' at  zero  when  a negative  value  (compresson  load)  would  normally 
J 

occur.  Reloading  follows  C-B-A.  The  general  load-deflection  curve  for  a 
compression-only  element  would  look  like  Figure  1-8  except  the  non-zero  loading 
would  be  in  the  lower  left  quadrant.  The  input  data  specifies  whether  the  beam 
is  a tension-only  or  compression-only  element,  if  it  has  been  identified  as  an 
unsymmetrical  element. 


In  order  to  model  a general  unsymmetrical  element  with  different  non-zero 
tension  and  compression  properties,  two  separate  beams  connecting  the  same  end 
points  are  used  in  the  model.  One  beam  is  defined  as  a tension-only  element, 
and  the  other  as  a compression  only  element. 


1.3. 5. 3.6  Damping  Forces  and  Moments 

In  addition  to  the  nonlinear  strain  forces  and  moments  derived  in  the  pre- 
ceding sections,  program  KRASH  also  provides  for  the  calculation  of  damping 
forces  and  moments  for  the  internal  beams.  In  order  to  retain  the  balance 
between  force  and  moment  for  the  coupled  degrees  of  freedom  inherent  in  the 
strain  formulation,  the  damping  forces  are  calculated  using  the  linear  stiffness 
matrix.  Naturally,  however,  these  forces  are  in  phase  with  the  velocities 
rather  than  the  displacements.  The  resulting  model  is  of  the  following  form. 


1fd[  *-|i[K]  jvel| 


(1-^48) 
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The  linear  stiffness  matrix  is  multiplied  by  the  beam  end  point  velocity 
vector.  For  pure  sinusoidal  motion,  this  formulation  is  the  same  as  conventioneU. 
structural  type  damping,  in  which  the  deunping  forces  are  proportional  in  magnitude 
to  the  strain  forces  but  in  phase  with  the  velocity.  In  such  a system,  w would 
be  the  response  frequency  and  ^ the  damping  ratio. 

For  the  actual  non-sinusoidal  motion  of  program  KRASH,  the  deunping  force 

model  could  be  considered  as  a pseudo-structural  damping  model.  In  this  case, 

the  w in  the  above  equation  is  taken  as  the  natiiral  frequency  of  the  isolated 

system  consisting  of  beam  ij  and  the  end  masses  m.  and  m, . The  actual  computa- 

^ J 

tional  sequence  used  in  program  KRASH  is  the  following. 

The  six  damping  forces  and  moments  at  beam  end  JN  are  calculated  using 
equations  having  the  same  form  as  (1-32),  (l-Ul)  and  (1-^*2),  depending  upon 
the  beam  end  conditions.  The  deflections  d^  - d^  in  those  eqviations  are  re- 
placed with  the  beam  velocities  vel^  - vel^  which  are  calculated  as  follows 

vx. 

1 

vy. 

VZ  . 

1 

v<t>. 

vSi 

v<|^ 

vG. 

1 

VlVj 

where  the  terms  on  the  right  are  given  by  Equations  (l-20)  and  (l-22).  Simi- 
larly, the  terms  dG^ , dGj,  d4>j,  d4>j  appearing  in  (l-Ul)  and  (1-U2)  for  fixed- 
pinned  end  conditions,  are  replaced  by  vG^,  vGj  , and  wi'j , respectively. 

The  resulting  damping  forces  at  JN,  in  beam  axes,  are  denoted  by  |fDj|,  and  are 
shown  here  for  reference. 


(1-149) 
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MD4ij  = 0 


(Pinned  at  j about  y axis) 

(l-51a) 


k.  k 

FDz.  = vel^  + v6. 


^ ^ ^ ) (Pinned  at  J about  z axis) 


MDGj  = 0 


(l-51b) 


FDy,  = 


MD4<  = 


FDz,  = 


MDe,  = 


J 


(Pinned  at  i about  y axis) 

(l-52a) 


Pinned  at  i about  z axis ) 
(l-52b) 


The  result  of  Equations  (1-50),  (l-5l),  and  (1-52)  is  a 6 element  force/ 
moment  vector  |fDj|.  Each  of  these  elements  is  next  multiplied  by  a damping 
coefficient  to  obtain  the  actual  deunping  force/moment  at  the  J end  of  the  beam. 


1-39 


(1-53) 


The  damping  coefficients  are  calculated  as 


=K-2V"k  k = l.  ....6 


(1-5M 


where  is  the  damping  ratio  for  beam  ij  which  is  an  input  constant.  The 
natural  frequency  for  each  of  the  6 beam  directions  is  calculated  aL 


V m^ 


mj) 


k = 1,  2,  3 


(I  + I ) 

=/  Ik  Jk  1,  = 1 o o 

k+3^/  r-^ ^ 1,  2.  3 

ik  Jk 


(l-55a) 


(l-55b) 


Equations  (1-55)  are  the  natural  frequency  eqviations  for  a beam  with  con- 
centrated inertias  at  each  end.  The  6 Kj^'s  used  in  (1-55)  are  Just  the  diagonal 
elements  of  the  stiffness  matrix  for  beam  iJ. 


(1-56) 
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These  terms  were  defined  in  Equations  (1-28),  (l-29c)  and  (l-30c).  The 
inertia  terms  and  k = 1,  2,  3,  are  the  inertias  of  mass  i and  J,  re- 

solved into  beam  ij  axes.  These  are  given  by  the  diagonal  elements  of  the 
following  matrices. 


= [aIJTAI]  [i^]  [aiJTAi]^  (l-5Ta) 

mass 
axes 

= [AIJTAJ]  [ij]  [aiJTAj]^  (l-57b) 

mass 
axes 

All  quantities  in  equations  (1-5^)  through  (1-57)  are  constants  which  do 
not  vary  with  time.  Therefore,  these  items  are  calculated  Just  once  prior  to 
the  main  time  loop  in  program  KRASH.  Equation  (l-53)  for  the  damping  force 
vector  is,  of  course,  within  the  time  loop.  In  Equation  (1-57),  the  time  zero 
values  of  [AIJTAI]  and  [AIJTAJJ  are  used,  since  these  are  functions  of  [A^^j], 

[A^]  and  [Aj].  Although  all  three  of  these  component  matrices  are  time-varying, 
their  time  zero  values  correspond  to  an  undistorted  vehicle  and  yield  the 
proper  transformation  from  mass  to  beam  axes  for  the  undeformed  structure. 

1.3. 5.3. 7 Total  Internal  Beam  Forces  and  Moments 

Equations  (l-l*6a)  and  (1-53)  yield  the  strain  and  damping  forces  and 
moments  at  JN.  These  are  denoted  jF  and  Jfd, respectively.  The  cor- 
responding  quantities  at  the  iM  end  of  the  beam  are  calculated  using  the  static 
balance  Equations  (l-32b).  These  are  denoted  { | *ind  These  four 

independent  sets  of  forces /moments  are  used  to  calculate  the  internal  beam 
strain  and  damping  energies.  These  calculations  are  described  in  Section  1.3.16. 
After  the  energy  calculations,  we  no  longer  need  the  strain  and  damping  forces 
separately.  Therefore,  the  two  are  now  added  together  to  form  the  total  inter- 
nal beam  load. 


iJ  axes 
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|FSi|=  |F^';  + jFD^'{ 


(l-58a) 


jFSj}  = jFj'l  + jFDj’  { 


(l-58b) 


These  forces  are  in  beam  axes  with  positive  directions  indicated  in 
Figiire  1-4,  which  shows  the  forces  acting  on  beam  ij.  To  get  the  forces  acting 


on  masses  m.  and  m , we  must  reverse  the  signs  and  resolve  the  vectors  from 
^ J 


beam  axes  to  mass  axes. 


jFS^'l  = - [AUTAl]"^  jFS^} 


(l-59a) 


|FSj’^  = - [AUTAJ]"^  |FSj} 


(l-59b) 


where  the  transformation  matrices  are  defined  by  Equation  (l-l8).  These  forces 


act  at  beam  ends  iM  and  JN  in  mass  axes. 


The  moments  acting  on  the  masses  must  now  be  corrected  to  account  for  the 
offset  of  the  node  points  from  the  masses.  Figure  1-2  shows  the  offset,  where 
the  vector  jrx,  ry,  rz[  contains  the  offset  components  in  mass  fixed  axes. 

These  are  not  time-varying  quantities,  since  each  node  point  is  rigidly  attached 
to  a defined  mass  point.  The  moments  at  the  mass  points  are  given  by 


K]  h: 


(l-60a) 


s r‘6 


(l-59a) 


(l-60b) 


where  the  first  vectors  on  the  right  hand  side  are  the  three  moments  from  Equa- 
tion (1-59).  These  moments  act  at  the  node  points  rather  than  at  the  masses. 
The  secohd  vector  contains  the  three  forces  from  (1-59)  which  act  at  the  node 
points.  The  transformation  matrices  are  of  the  following  skew-symmetric  form: 


-rz  ry 

0 -rx 

rx  0 


iM 


(l-6la) 


These  are  not  time- varying  matrices. 


(l-6lb) 


Equations  (l-59)  and  (l-6o)  give  the  current  internal  beam  forces  and 
moments  acting  on  masses  m^  and  m^ , due  to  beam  ij , in  mass-fixed  axes.  To  ob- 
tain the  total  internal  forces  and  moments  acting  on  each  mass , the  contribu- 
tions from  each  internal  beam  connecting  to  a given  mass  must  be  accounted  for. 
In  program  KRASH,  this  is  done  within  the  beam  iJ  loop  where  | andjFSj'  | 

are  calculated.  These  are  merely  added  to  the  total  internal  forces  emd  moments 
for  mass  i and  mass  J . 
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(l-62b) 


Equations  (1-62)  should  be  interpreted  as  follows.  At  each  time  cut,  prior 
to  entering  the  beam  loop  which  calculates  the  beam  forces  for  all  the  ij  beajns 
in  the  model,  (Xj^  ...,  are  all  set  to  zero  for  all  masses  i.  Equation 

(1-62)  is  then  employed  for  each  beam  in  tvirn,  where  andjPSj'f  are  added 

to  the  total  internal  beam  forces  acting  on  the  proper  masses  i and  j , respec- 
tively. Once  all  the  beams  have  been  calculated  in  the  beam  loop,  ..., 

Nj^)  will  contain  the  contributions  from  all  beams  connecting  to  mass  m^,  for 
all  the  masses  in  the  model. 


1.3.5***  External  Forces  and  Moments 

The  crash  forces  result  from  the  compression  of  external  springs  radiating 
from  selected  mass  or  node  points,  when  these  springs  contact  either  the  ground 
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plane  (z  = 0)  or  the  slope  = o).  The  following  sketch  shows  the  three 

spring  ladlating  outward  from  mass  i; 


The  free  (uncompressed)  lengths  of  the  three  springs  are  denoted 

i and  % ..  These  springs  lie  along  the  body  ajces  for  the  ith  mass.  JJ, 
yi  zi 

and  are  input  either  positive  or  negative,  depending  upon  whether 
the  spring  radiates  outward  along  the  positive  or  negative  body  axis.  Those 
shown  above  are  all  positive  springs.  points 

of  the  fully  extended  springs.  Only  data  for  those  springs  desired  in  the 
analysis  need  be  input;  if  no  spring  is  input,  its  crash  forces  are  set  to 
zero.  Typicallj',  only  those  masses  at  the  extremities  of  the  vehicle  likely 
to  contact  the  ground  would  have  crash  springs.  The  intent  of  these  springs 
is  to  represent  crushable  external  structure  such  as  found  on  the  lower  fuse- 
lage and  nose  sections. 

The  first  step  in  the  calciilation  is  to  determine  which  of  the  C . , C . , 

xi  yi 

and  spring  end  points  are  below  the  ground  or  slope.  Only  these  springs 
are  further  Malyzed  to  determine  their  crash  forces.  For  these  springs,  the 
spring  compression  along  the  spring  axis  is  determine.’,  as  well  as  the  com- 
pression velocity  and  the  ground  contact  point  velocity.  Then  the  spring 
force  along  the  spring  axis  is  calculated,  from  which  the  vertical  load  in 
ground  axes  at  the  ground  contact  point  is  obtained.  From  this  vertical  load 
and  an  input  friction  coefficient,  the  ground  drag  load  is  computed.  The  direc- 
tion of  this  load  is  opposite  the  ground  contact  point  velocity  vector. 


Finally,  the  three  forces  at  the  ground  contact  point  are  rotated  into  ith 
mass  body  axes,  and  the  resulting  forces  and  moments  at  mass  i are  calculated. 


1.3. 5.^.1  Spring  Compression 

Figure  1-9  represents  a typical  spring  contacting  the  ground  at  an  oblique 
angle.  This  figure  illustrates  the  most  general  case  where  the  spring  is 
attached  to  a node  point  rather  than  a mass  point.  External  springs  can  be 
attached  either  directly  to  mass  points  or  to  node  points.  Point  is  the 
ground  contact  point.  Treating  the  fully  extended  spring  as  a vector,  that 
portion  of  the  vector  below  the  ground  represents  the  spring  compression. 

This  is  the  distance  from  point  C!  to  point  C^.  The  coordinates  of  point 
in  ground  axes,  are  given  simply  by 


xc  . xc  . xc  . 
XI  yi  zi 
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xl  yl  zl 
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rx  rx  rx 


rz  rz  rz 


The  above  equation  gives  the  coordinates  of  all  three  for  mass  i; 

these  points  are  denoted  C . , C . , and  C . . Figure  1-9  shows  only  one  of  the 

XX  Zl 

three  springs,  so  the  nomenclature  in  Figure  1-9  is  not  subscripted  to  indi- 
cate whether  an  x,  y or  z spring  is  shown.  The  third  term  in  the  above 
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The  dxc^,  dyc^,  and  dzc^  in  Figure  1-9  can  be  determined  from  proportion- 
ality considerations  as 
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and  similarly  for  the  y and  z springs, 
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From  these  ground  axes  components  of  the  vector  from  to  , the 


magnitude  of  this  vector  (equal  to  the  spring  compression)  is  given  simply  by 
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In  order  to  simplify  the  notation,  the  following  vectors  and  matrices  are 
iefinei  'i  = 1,  indicates  the  mass;  J = 1,2,3  indicates  the  direction; 
k = 1.2,3  indicates  which  spring): 


The  above  equations  are  Appropriate  if  the  spring  intersects  the  horizontal 
ground  plane.  If  the  spring/intersects  the  sloping  ground  plane,  the  following 
equations  are  utilized: 
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vcstj3j  - (sin  e)  ♦ (co.  6) 


(1-66) 


I.3.5.U.2  Ground  Contact  Point  Velocity^ 

The  next  computation  required  is  to  determine  the  spring  contact  point 
velocity  in  the  ground  plane  (either  the  horizontal  ground  plane  or  the  slope 
plane,  whichever  the  spring  is  contacting).  This  is  necessary  to  detemine 
the  line  of  action  of  the  ground  plane  drag  load,  which  is  opposite  to  the 
contact  point  velocity  vector.  The  velocity  of  point  C with  respect  to 
ground  is  given  in  vector  notation  as 


0 -*C'  m.  — -C  0 -*C' 

V = ^ V + V 


(1-87) 


where 


0 -■  » C ' 

^ = velocity  of  C*  with  respect  to  ground 


m.— C = velocity  of  C’  with  respect  to  m. 

1 V 

0 - ■»  c ' = velocity  of  point  fixed  in  m. , instantaneously  coinciding  with 

V point  C,  with  respect  to  ground 


This  latter  term  is  given  by 
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velocity  of  with  respect  to  ground 


w = rotational  velocity  of  with  respect  to  ground 
^IkM  * vector  from  to  spring/ grovind  contact  point  C* 

Figure  1-10  illustrates  the  geometry  involved.  The  vector  R.,..  is  given  by 
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The  components  of  (1-8T)  and  (1-88)  are  given  by 
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Figure  1-10.  External  Spring  Geometry 


Combining  the  above  equations,  the  ground  contact  point  velocity  is  given 
by 


- 0 -r^  q^- 
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(1-93) 


where  the  three  options  apply  to  k = 1,  2 or  3 In  the  order  shown.  In  other 
words.  Equation  (1-93)  represents  three  equations,  applicable  to  either  an  x, 
yorz(k=l,  2or3)  direction  external  spring.  The  time  varying  spring 
length  l^j^  is  Equation  (1-93)  is  calculated  as 


‘ikM  = 


(I-9I1) 


This  form  is  used  so  that  has  the  same  sign  as  the  input  which 

as  mentioned  earlier  can  be  negative  to  allow  springs  in  the  negative  direc- 
tion of  the  body  axes.  These  spring  lengths  are  shown  in  Figure  1-10;  the  points 
C in  this  figure  are  the  ground  contact  points.  If  an  external  spring  is 
directly  connected  to  a mass  m^^  instead  of  a node  point,  then  (rx,  ry,  rz)  are 
all  zero. 
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If  the  spring  is  contacting  the  slope  rather  than  the  horizontal  ground, 
the  above  velocities  are  transformed  into  slope  axes  as  follows: 


(1-95) 


1.3.5*^-3  External  Spring  Axial  Force 

This  section  discusses  the  calculation  of  the  externeLL  spring  force  along 
the  spring  axis.  The  axial  force  in  the  spring,  denoted  is  input  as 

a tabular  function  of  the  spring  compression  As  with  the  internsJ.  forces, 

unloading  along  em  elastic  line  is  included.  However,  unlike  the  internal 
forces,  an  extension  load  is  not  allowed;  positive  is  a conjjressive 

load  acting  at  the  ground  contact  point  along  the  axis  of  spring  ikM. 

Figiire  1-11  illustrates  the  shape  of  the  input  external  spring  load- 
deflection  curve.  SI,  SA,  SB,  SF,  FSPOI,  FSPOF  and  Ke  are  all  input  constants. 


SPRING  AXIAL 
FORCE 


FSPOF,, 


FSPOI 


ikM  B. 


V l»L-/ 

b s,  Q 


Sa  H 


SPRING  AXIAL  COMPRESSION 

Figure  1-11.  External  Spring  Load-Stroke  Curve 
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Loading  proceeds  eLLong  0-A-B-C-D-E.  At  any  spring  compression  less  than 

SA,  unloading  occurs  along  a line  having  the  same  slope  as  0-A.  At  any  spring 
compression  between  SA  and  SF,  unloading  occurs  along  a line  parallel  to  B-C. 
(If  an  extension  of  this  unloading  line  would  cross  above  0-A-B,  resulting  in 
more  energy  released  dui'ing  unloading  than  input  during  loading,  then  the 
xmloading  curve  is  parallel  to  0-A  rather  than  B-C. ) Once  SF  is  reached,  all 
subsequent  loading  and  unloading  follows  line  D-E  emd  its  extension  in  both 
directions.  In  all  cases,  unloading  does  not  extend  beyond  FSPO  = 0.  For 
example,  if  unloading  begins  at  F,  it  proceeds  along  F-G  (parallel  to  0-A), 
and  then  holds  at  FSPO  = 0 from  G to  H.  Subsequent  reloading  takes  place  eulong 
H-G-F-B-C-D-E.  There  is  no  rupture  of  external  spring  elements. 

A first-order  approximation  of  the  effect  of  ground  flexibility  is  in- 
cluded in  the  following  manner.  A ground  flexibility  constant  is  input 

for  each  external  spring.  The  deflection  coordinates  in  Figure  1-11  are  modi- 
fied to  include  this  flexibility. 


S^^^ikM  = S^ikM  ^ FSPOI^j^^  * fg^^ 


" ®^ikM  * ^^°^ikM  * ^®ikM 


S^^^ikM  = S®ikM  ^ ^P°^ikM  * "SikM 


®^^ikM  ■ ^^ikM  * ^®^°^ikM  * ^®ikM 
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This  has  the  effect  of  softening  the  input  load-deflection  curve  representing 
the  airplane  structure  to  indlude  the  flexibility  of  the  ground  in  series  with 
the  structure.  The  bottoming  spring  constant  is  modified  as  follows: 


The  ground  deflection  is  calculated  for  each  spring  at  each  time  cut  as 


DELGi^  = FSPO^j^  * fg^j^ 


(1-98) 


Furthermore,  the  ground  is  never  allowed  to  rebound;  the  ground  deflection 
always  increases  or  stays  constant.  This  is  accomplished  by  setting  DELG^j^ 
at  each  time  step  to  the  maximum  of  the  current  value  and  the  previous  maxi- 
mum. The  ground  deflection  must  also  be  accounted  for  in  the  calculation  of 
the  deflected  spring  length  from  Equation  (1-9^). 


'ikM  ' 


(^ikM) 


3.,  ,,  + DELG.,  „ 
ikM  ikM 
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In  effect,  the  spring  compression  is  reduced  by  the  amount  that  is 

directly  attributable  to  ground  deflection,  since  the  spring/ groiond  interface 
is  actually  lower  than  would  be  calculated  for  a rigid  ground. 

The  flow  diagram  for  the  calcxilation  of  the  spring-force  FSPO^j^  is  shown 
in  Figure  1-12.  The  calculations  are  performed  within  nested  i and  k loops. 
Immediately  prior  to  the  calculations  shown  in  Figure  1-12,  there  is  a test 
to  see  if  an  external  spring  exists  for  the  Ik  pair  being  considered.  If  a 
spring  exists,  it  is  checked  to  see  if  it  is  in  contact  with  either  the 
horizontal  ground  or  the  slope.  If  both  these  tests  are  satisfied,  the 
algorithm  shown  in  Figure  1-12  is  entered,  otherwise  this  block  of  coding  is 
skipped.  In  Figure  1-12,  SK  is  the  spring  compression  s^j^. 


The  test  on  SDOT  determines  whether  loading  or  unloading  is  occurring. 
However,  if  SDOT  is  greater  than  zero,  this  may  be  either  loading  (as  from 
A to  B in  Figure  1-11 ) or  reloading  (as  from  G to  F in  Figure  1-11 ).  In 
order  to  distinguish  between  these  situations,  a test  is  performed  to  see  if 
SK  is  greater  than  SBAB.  SBAR  is  initially  set  to  zero,  and  during  loading 
it  is  always  equal  to  the  spring  compression  SK  from  the  previous  time  step. 
Therefore  during  loading  this  test  is  passed,  and  the  path  is  followed  wherein 
FSPO  is  calculated  directly  from  the  input  data. 

In  this  coding  path,  IBS  is  the  bottoming  spring  flag.  This  is  initially 
set  to  zero,  and  remains  zero  until  the  spring  compresses  beyond  point  D in 
Figure  1-11.  Once  the  bottoming  spring  is  reached  (SK  ^ SFFL) , IBS  is  set  to  1 
and  all  subsequent  calculations  of  FSPO  use  the  equation  representing  line 
D-E  and  its  extension  in  Figure  1-11.  In  the  basic  loading  path,  FSPBAR  and 
SBAR  are  set  to  the  current  values  of  spring  load  and  deflection. 

Once  unloading  is  initiated,  SDOT  will  be  less  than  zero.  Then  the 
counter  NN  is  tested;  if  this  is  zero  then  this  is  the  first  time  step  of  the 
current  unloading.  The  bottoming  spring  flag  IBS  is  checked,  SK  is  compared 
with  SA,  the  appropriate  equation  is  used  to  calculate  the  xinloading 
stiffness  Kl)N,  and  NN  is  set  to  1.  FSPO  is  then  calculated  using  the  unloading 
stiffness,  SBAR  and  FSPBAR,  which  are  the  compression  and  load  when  unload- 
ing began.  For  example,  SBAR  and  FSPBAR  are  the  coordinates  of  point  F in 
Figure  1-11  where  unloading  began.  Subsequent  unloading  paths  bypass  the 
calculation  of  the  unloading  stiffness,  since  NN  is  now  1. 

Once  SDOT  again  becomes  positive,  SK  is  initially  less  than  SBAR  emd  re- 
loading occurs,  which  uses  the  same  equation  as  unloading.  Finally,  after 
sufficient  reloading,  SK  will  be  greater  than  SBAR,  and  the  basic  loading  path 
is  again  followed.  In  this  path,  NN  is  reset  to  zero  and  SBAR  and  FSPBAR  are 
continually  updated.  Finally,  regardless  of  the  method  of  calc\ilating  FSPO, 
if  it  is  calciilated  as  less  than  zero  it  is  reset  to  zero. 

1.3.5.^*^  Crash  Forces  at  Ground  Contact  Point 

The  three  components  of  the  grovind  interaction  force  at  the  ground  con- 
tact point  are  denoted  by  XVOC  , where,  as  before,  i refers  to  the  ith 

1 JX 
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mass,  J to  the  force  direction,  and  k to  the  kth  spring  for  mass  i.  First 
the  component  of  the  axial  spring  force  perpendicular  to  the  ground  is 
computed.  This  is  positive  upward  acting  on  the  spring.  This  is  given 
simply  by  resolving  into  ground  axes  and  retaining  only  the  vertical 


ccmoonent. 


XVOC.„  = A.  „FSPO,.  SIGN 
i3k  i3k  ik 


('iJ 


k = 1,2,3 
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since  FSPO^j^  acts  along  the  kth  body  axis  for  mass  i.  The  SIGN  term  is 


necessary  to  retain  the  convention  of  XVOC.  acting  positive  upward  with 

1 


negative  springs.  The  ground  plane  components  of  the  axial  force  FSPO  are 

1.K 


ignored,  since  these  values  are  next  computed  on  the  assumption  that  they 
are  functions  of  the  ground-spring  friction  coefficient.  The  ground  plane 
velocity  vector,  from  Equation  (l-93),  has  the  magnitude 


ik 


= [(vcp^ 


Ik)^  + (^<^Pi2k^^] 


1/2 
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This  is  shown  in  the  following  sketch,  a plan  view  looking  downward  at 


ground  contact  point  C : 

ik 


r 


The  nagnltude  of  the  resultant  ground  plane  load  is  given  by  (mu.,XVOC. 

llv 


^piowiit)*  i*®*»  input  friction  coefficient  times  the  vertical  load 
perpendicular  to  the  ground  plane,  plus  a constant  plowing  force,  ^plowik 
and  are  input  constants.  Fpiowij^  constant  plowing  force  that  acts 

from  t = 0 to  t = PLOWT,  an  input  constant,  and  then  drops  to  0 for  the 
remainder  of  the  run.  and  ^piow^j^  together  allow  the  program  user  to 

calculate  a drag  load  that  is  either  dependent  on  the  vertical  load  or  is  a 
constant.  The  direction  of  the  total  ground  plane  load  is  opposite  to  the 
contact  point  velocity  vector  as  shown  in  the  above  sketch.  Therefore, 


the  drag  and  side  loads  at  point  in  ground  eixes,  are  given  by 


= [■ 


XVOC...  = mu,,XVOC.,,,  + F , 
ijk  ik  i3k  plow 


ik 
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These  are  illustrated  in  the  above  sketch.  These  forces  are  now  resolved  back 
into  ith  mass  body  sixes;  these  resolved  forces  are  denoted  by  FSP^ 

We  have 
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The  minus  sign  is  used  so  that  these  forces  have  the  same  sense  as  the 


basic  body  axis  convention,  since  XVOC.  „ were  positive  aft,  left  and  upward. 

X jK 


These  forces  are  illustrated  in  Figure  1-13. 


If  the  spring  contacts  the  inclined  slope  rather  than  the  horizonteLL 


ground,  €?.l  three  ground  axis  components  of  the  axial  spring  force  FSPO 


ik 


are  obtained  first,  and  then  resolved  into  slope  axes,  retaining  only  the 
component  normal  to  the  slope.  For  this  situation. 


J 
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rcea  and  Momenta 


(1-lOU) 


from  which  the  con5)onent  normal  to  the  slope  is  obtained  as 
XVOCSL...  = (sin  6)  XVOC.,.  + (cos  B)  XVOC 

ixK  X^K 


(1-105) 


where  8 is  the  slope  angle  measured  from  horizontal. 

The  contact  point  velocity  vector  in  the  inclined  plane,  from  Equa- 
tion (1-95),  has  the  magnitude 
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The  drag  and  side  loads  in  the  plane  of  the  inclined  slope  are  as 
follows . 
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muik  XV0CSL^3j^ 
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Finally,  the  three  contact  point  forces  given  by  Equations  (1-105)  and  (1-107) 
are  resolved  into  horizontal  ground  axes  and  then  into  ith  mass  axes  to  form 
the  FSP^jj^  illustrated  in  Figure  1-13. 
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k = 1,2,3  (1-108) 


Now  that  the  spring  contact  forces  are  in  ith  mass  axes,  there  is  no 
further  distinction  between  forces  resulting  from  contact  with  the  horizontal 
ground  (Equation  (1-103))  and  those  resulting  from  contact  with  the  slope 
(Equation  (I-IO8)). 

1.3. 5. **.5  Crash  Forces  and  Moments  at  Point  mj 

Figure  1-13  shows  the  crash  forces  and  moments  acting  at  m^ , which  are 
Just  the  resultant  of  the  FSP...  also  shown  in  Figure  1-13;  i.e.,  the  FSP. 

1>JK.  IjK 

are  just  moved  to  point  m^  and  compensating  moments  at  m^  added. 
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The  are  from  Equation  (1-99);  recall  that  they  may  be  plus  or  minus  to 
allow  for  springs  in  the  negative  direction  of  the  mass  body  axes.  Figure  1-13 
shows  all  positive. 

Equations  (1-109)  yield  only  the  contribution  of  springs  attached  to  node 
point  iM  to  the  total  external  spring  forces  and  moments  acting  on  mass  m^. 

To  achieve  the  desired  total  forces  and  moments  on  m^,  the  contributions  from 
all  node  points  attached  to  m^  and  having  external  springs  must  be  summed. 


(1-110) 


where  NNP^  is  the  number  of  node  points  attached  to  mass  m^  and  having  external 
springs . 

1.3.6  Total  Forces  and  Moments 

The  total  forces  and  moments  at  m^^ , in  ith  body  axes,  are  used  to  drive 
the  six  rigid-body  degrees  of  freedom  of  lumped  mass  m^.  These  total  forces 
and  moments  are  Just  the  sum  of  the  components,  i.e.,  the  gravity,  aerodynamic, 
internal  beam,  and  external  spring  forces  and  moments.  The  total  loads  are 
given  by: 


J 


There  are  no  gravity  and  aerodynamic  moments.  The  terms  of  Equation 
( 1-111)  come  from  Equations  (1-5) » (1-8),  (1-62)  and  (l-llO). 

1.3.7  Rigid-Body  Equations  of  Motion 

Having  obtained  all  the  forces  and  moments  acting  on  each  mass  m^ , the 
rigid-body  equations  of  motion  for  mass  m^  can  be  written.  Euler's  equa- 
tions of  motion  are  used.  These  are  derived  in  Reference  1;  the  derivation 
is  not  repeated  here.  The  body  axis  components  of  the  absolute  (relative  to 
ground)  translational  velocity  of  mass  m^  are  denoted  by  u^,  v^,  and  w^.  The 
corresponding  rotational  velocities  are  designated  p^ , q^ , and  r^ . These 
velocities  are  shown  in  the  following  sketch: 

I 

body  axes  fixed  in 
i 


Z-iler's  equations  of  mction  are  six  equations  for  the  time  derivatives 
of  the  six  velocities,  in  terms  of  the  six  total  forces  and  moments  from  the 
preceding  section,  the  velocities  themselves,  and  the  inertia  properties  of 
mass  m^.  These  equations  are  (from  Reference  l): 
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X. , Y.,  Z.,  L.,  M. , and  N.  are  the  total  forces  and  moments  from  Equation 
i i i i i 1 

( 1-111 ).  I .,  I and  I . are  the  moments  of  inertia  of  mass  m.  about  its 
XI  yi  zi  1 

body  fixed  axes.  I .,  I .,  and  I . are  the  products  of  inertia,  e.g. , 

''  xyi  yzi’  xzi  ^ 

I . = / x.y.dm. , etc.  H . , H . , and  H . are  the  angular  momenta  of 
xyi  11  1 exi  eyi  ezi 

masses  m^,  due  to  rotation  of  internal  masses  within  m^ , such  as  propellers, 
rotors,  and  engine  turbines. 

The  above  equations  of  motion  are  integrated  numerically  to  yield  u^, 

^i*  '^i’  ^i’  ^i‘  translational  velocities  are  obtained  in  ground 

axes  by  a simple  transformation. 


li  } - [Ajl  ^ Vj 


(1-lli*) 


Similarly,  the  angular  velocities  p^,  q^,  and  r^  are  the  bo(3y  axes  com- 
ponents of  the  instantaneous  absolute  rotation  rate  of  mass  ra^.  These 
variables  do  not  correspond  to  any  set  of  coordinates  which  specify  the 
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orientation  of  the  airplane.  Therefore,  in  order  to  solve  for  the  orienta- 
tion, it  is  necessary  to  transform  the  p',  q*,  r'  to  the  time  derivatives 
of  the  Euler  angles.  By  summing  the  orthogonal  projections  of 
*^i’  ®i*  '^i  fixed  axes,  the  following  kinematic  relations 

are  obtained  (Reference  1): 


+ (q^8in(t)^+rj^co8(J)^)  tan  0^ 

®i  “ ■ rj8in(()^  (1-115) 

<l'i  ® (q^8lni^j^+rj^cos(J>^)  8ec  0^^ 


Equations  (l-llU)  and  (1-115)  are  integrated  numerically  to  obtain  x^,  y^^, 
and  4>^,  0^,  4;^.  Also,  p^^,  q^,  and  r^  are  integrated  to  obtain  the  quan- 
tities Inp^,  Inq^,  and  inr^.  The  Incremental  changes  in  these  Integrals  are 
used  in  Equation  (l-12).  They  represent  the  incremental  rotations  of  mass  m^^ 
in  body  fixed  axes. 

It  should  be  noted  that  (li^,  v^,  w^)  from  Equation  (1-112)  are  not  the 
components  of  the  absolute  acceleration  of  m^^.  The  expressions  for  these  quan- 
tities are; 


(1-116) 


where  the  accelerations  are  in  g units.  These  are  used  only  for  output  data  in 
program  KRASH.  Con?)aring  Equations  (1-112)  and  (I-II6),  the  relationships 
between  the  absolute  accelerations  and  the  time  derivatives  of  the  body  eucis 
components  of  the  absolute  velocity  are  seen  to  be: 
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These  relationships  can  be  independently  derived  from  the  time  derivative 
of  the  absolute  velocity,  accounting  for  the  time  rate  of  change  of  the  unit 
vectors  fixed  in  the  body. 

1.3.8  Filtered  Accelerations 

Program  KRASH  has  been  extensively  used  in  conjunction  with  full  scale 
crash  tests , wherein  it  is  desired  to  establish  the  degree  of  correlation 
between  analytical  and  test  results.  The  most  readily  obtainable  test  results 
are  measured  accelerations  at  various  points  on  the  structure.  Since  crash 
test  accelerometer  outputs  generally  contain  a great  deal  of  very  high  fre- 
quency, high  amplitude  responses,  the  data  is  normally  filtered  prior  to  its 
usage  or  formal  presentation.  A wide  variety  of  analog  and  digital  filters 
may  be  employed,  but  in  general  the  cutoff  frequency  is  in  the  area  of  100  Hz 
to  300  Hz. 

In  order  to  facilitate  comparisons  of  KRASH  analytical  results  with  test 
results,  the  calculation  of  filtered  mass  accelerations  is  included  in  the 
program.  The  filter  grapioygd  .in  the  program  is  a simple  first  order  lag.  The 
Laplace  transform  for  such  a filter  is  given  by: 

1 

T(‘**)  = I ■ U)  (1-118) 

where 

j = 

^ = response  frequency,  rad/sec 

w ■ cutoff  frequency  of  filter,  rad/sec 
c 
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The  phase  and  amplitude  response  for  this  filter  is  shown  in  Figure  1-lU. 
Although  the  filters  used  for  test  data  reduction  are  generally  much  steeper, 
it  is  felt  that  the  computational  effort  (and  run  time)  required  for  a more 
complex  filter  is  not  appropriate  for  inclusion  in  the  direct  time-history 
solution  of  program  KRASH. 

The  first  order  filter  is  implemented  throvigh  the  following  equation  in 
the  time  domain: 


py  + y = X 


(1-119) 


where 


X = unfiltered  acceleration  (input)  time  history 

y = filtered  acceleration  (output)  time  history 

y = y(t)  (1-120) 

- i _ 1 

^ “ WC  ~ 2TTfc 


The  filter  cutoff  frequency  fg  is  an  input  peirameter. 
Equation  (1-119)  is  solved  for  y: 


y = ^ (1-121) 

This  derivative  is  then  integrated  directly  along  with  the  other  integrations 
performed  in  the  program  to  yield  y(t)  which  is  the  desired  filtered  acceleration 
This  procedure  is  employed  for  all  the  linear  accelerations  for  all  the  masses 
in  the  model.  The  unfiltered  accelerations  used  for  the  input  to  (1-121)  are 
those  from  Equation  (l-ll6),  which  are  in  g units.  Therefore,  the  filtered 
accelerations  are  also  in  g units. 
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1.3.9  Translational  Momentum  Equations 


r 


Once  the  velocities  of  all  the  lumped  masses  are  known  from  the  pre- 
ceeding  equations,  it  is  possible  to  calculate  the  vehicle  eg  velocity  based 
on  momentum  considerations.  This  is  done  only  for  the  purpose  of  providing 
additional  output  information  for  the  program  user;  the  eg  velocities  thus 
calculated  are  not  used  elsewhere  in  the  program  computations.  From  conserva- 
tion of  linear  momentum; 
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(1-122) 


where 


_AP| 

^ AP  I 

^AP 


X. 
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^i 


= ground  axis  components  of  overall  airplane  eg  velocity 


W 


w • 

TOT 
N 


= ground  axis  components  of  ith  mass  translational  velocity 

= weight  of  ith  lumped  mass 
= total  vehicle  weight 
= number  of  lumped  masses 


The  results  of  Equation  (1-122)  are  printed  at  every  print  time  interval. 
This  information  is  useful  for  giving  the  program  user  an  indication  of  the 
translational  velocity  state  of  the  vehicle  as  a whole,  which  is  very  diffi- 
cult to  deduce  from  inspection  of  the  velocities  of  all  the  lumped  masses . 

A similar  calculation  could  be  performed  for  the  angular  velocity  of  the 
vehicle.  However,  since  the  angular  inertia  properties  in  ground  axes  are 
dependent  on  mass  orientation,  the  moments  of  inertia  for  all  the  lumped 
masses  would  have  to  be  continually  calculated  in  ground  axes,  as  would  the 
angular  velocities  of  the  masses,  since  their  spatial  attitude  varies  with 
time.  Alternatively,  the  angular  velocities  and  inertias  could  all  be  cal- 
culated in  some  reference  frame  approximating  the  overall  vehicle  attitude. 

In  either  case,  it  is  felt  thac  the  complexity  of  the  necessary  calculations 
outweighs  any  potential  benefits  from  the  overall  airplane  angular  velocities 
that  could  be  derived.  Therefore,  average  airplane  angular  velocities  are 
not  calculated. 

1.3.10  Control  Volume  Mass  Penetration  Calcxxlations 

The  computer  program  includes  the  calculation  of  whether  or  not  any  of 
the  N lumped  masses  have  penetrated  into  an  input  control  volume.  The  purpose 
of  these  calculations  is  to  determine  if  a major  mass  item  (such  as  the 
engine)  has  moved  into  a position  where  it  threatens  the  vehicle's  passengers. 
Therefore,  the  control  volume  input  is  intended  to  define  the  volume  of  the 
vehicle  in  which  human  occupants  are  present.  The  rectangular  control  volume 
is  located  with  respect  to  one  of  the  lumped  masses,  which  is  specified  in  the 
input.  Figure  1-15  shows  a typical  rectangular  control  volume  defined  with 
respect  to  mass  m^. 

The  six  walls  of  the  control  volume  are  always  perpendicular  to  the  three 
body  axes  fixed  in  mass  mp.  Thus,  the  spatial  orientation  of  the  control  vol- 
ume varies  as  mass  m rotates,  x , x , y , y , z , and  2 , all  positive, 
are  the  input  constants  defining  the  distances  from  m to  the  six  walls, 
measured  in  the  positive  and  negative  directions  along  the  body  fixed  axes 


of  m . 


CCHTROIi  KASS  n 


CCRTROL  VOUME 


Figure  1-15.  Control  Volume  for  Mass  Penetration 


Mass  penetration  is  determined  by  merely  examining  the  pth  body  axes  com- 
ponents of  a vector  from  control  mass  m^  to  mass  m^,  i = 1,2,,,,,N  (i  ^ p>.. 
These  components  are  given  simply  by 


■rfhere  [Ap3  is  the  transformation  matrix  [a^]  for  i * p.  These  vector  com- 
-ner.ti  are  then  tested  against  the  control  volume  dimensions.  If  all  three 


f-.f'onents  ile  within  the  -i  rrespondln#  pairs  -f  control  w . *e  walls 


AO-A055  696 
UNCLASSIFIED 


LOCKHEEO-CALIFORNIA  CO  BURBANK  F/6  1/2 

6ENERAL  AVIATION  AIRPLANE  STRUCTURAL  CRASHWORTHINESS  USER'S  MAN~ETC(U) 
FEB  76  M A 6AM0N  OOT-FA75WA-3707 

LR-26307  FAA/R0-77/169-V0L-1  NL 


-■S . 

t 

.4 

f 1.  4 . ' 

.- 

• 1 

- 

■ firr- 

P 

i'”  ' 

anTVZT 

•luy 

W' 

bifzt; 

, ..  ■ 

1 

mmmm  m > 

i 

■'u 

• 

M 

'iSk 

.■1‘ 

|Ayi. 

■ ■ 1 1 a-  • • 

it! 

k 

''■  ^ ' I*  •• 

ilijiiuiilh. 

'*  * ” 

1 ^ 

k 

IIM 

aWs 

■Lki 

1 i 
• ■ 

»±i 

»M  •.  ' ■ ■ I 

.1  ■ * ' 

!•••  ■ h h 

' * 

4i« . * . ? 

. 1. 

^Ukl 

■B 

Mijm! 

1 1 

liLilK: 

.fcj^ 

ilUl: 

£^1 

V4ju.i 

then  the  ith  point  has  penetrated  Inside  the  control  volume.  Whenever  a mass 
penetration  occurs,  the  program  prints  out  the  time  and  mass  that  penetrated. 
These  data  are  repeated  in  a summary  table  at  the  end  of  the  run,  showing  all 
the  mass  penetrations  for  that  run. 


1.3.11  Volume  Cheinge  Calculations 

In  order  to  determine  the  potential  for  occupant  injury  in  a crash,  the 
program  includes  the  capability  to  calculate  the  change  in  volume  of  a region 
defined  by  eight  (8)  corner  point  masses.  In  practice,  the  volume  specified 

would  represent  the  region  in  vdiich  human  occupants  are  located.  Figure  I-I6 
illustrates  a typical  volume  specified  by  defining  the  mass  numbers  at  the 
corner  points  1-8. 


2 


Figure  I-16.  Volume  Change  Geometry 


The  program  calculates  the  x,  y,  z ground  coordinates  of  the  circled 
points  @ through  © as  the  average  of  the  coordinates  of  the  appropriate 
I*  mass  points.  For  example,  to  calculate  the  coordinate  of  ©,  the 
coordinates  of  masses  located  at  1,  2,  3 and  1*  are  used.  When  the  ground 
coordinates  of  (l)  through  © are  known,  the  lengths  of  lines 
(D  - © and  © - © are  computed.  These  are  then  multiplied  to  obtain 
the  volume.  Both  the  volxime  and  the  lengths  of  lines  ©-©,©-©, 
and  © - © are  printed  at  each  print  time  interval.  The  program  allows  for 
up  to  6 separate  volumes  to  be  specified,  so  that  different  occupant  regions 
can  be  separately  einalyzed. 

While  Figure  l-l6  implies  that  8 masses  must  be  located  in  the  model  in 
a perfectly  rectang\ilar  pattern,  sloping  planes  are  acceptable  since  the 
volinne  is  calculated  using  coordinates  of  the  center  points  of  each  side 
plane.  Therefore,  in  any  three  dimensional  structural  model  there  should  be 
8 masses  located  sufficiently  close  to  the  corners  of  the  occupiable  volume, 
so  that  extra  mass  points  need  not  be  added  to  the  model  solely  for  volume 
change  calculations. 

1.3.12  Dynamic  Response  Index  (PRl) 

The  program  also  includes  the  capability  of  computing  Dynamic  Response 
Index  (DRI)  for  certain  selected  masses.  DRI  is  a nondimensionalized  measure 
of  the  compression  of  the  human  spinal  column,  and  statistical  data  is  avail- 
able relating  DRI  to  the  probability  of  spinal  injury.  References  2 and  3 
describe  DRI's  in  detail.  Basically,  all  that  is  involved  is  the  mathematical 
model  illustrated  in  Figure  1-lT- 

Mass  m^  is  the  occupant's  seat;  mass  m^  is  the  upper  torso  of  the 
occupant.  The  frequency  and  damping  ratio  given  are  intended  to  properly  model 
the  dynamic  response  to  the  upper  torso  to  the  seat's  excitation.  Mass  m^  is 
also  the  upper  torso,  but  the  stiffness  and  damping  are  slightly  different. 

This  is  done  to  obtain  a response  number  (DRI)  for  mass  m^  that  correlates 
veil  with  statistical  injury  data.  The  important  point  is  that  mass  m^  is 
driven  by  mass  m^^,  but  the  interconnecting  forces  only  drive  mass  m^,  not 
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Figure  1-17.  DRI  Model 


mass  m^.  Thus,  the  coupled  systems  of  m^  and  m2  determine  the  response  of 
m^,  and  this  response  drives  m^.  The  Dynamic  Response  Index  is  defined  in 
Figure  1-17  hy 


2 

0l>  6 

DRI  = -2-  (l-12h) 


Substituting  the  constants,  the  DRI  is: 

DRI  = 6.556  (1-125) 


If  the  damping  were  zero,  the  DRI  would  be  the  acceleration  of  mass  m^. 
With  damping,  it  is  a nondimensionalized  measure  of  spring  compression.  In 
the  computer  program,  6 is  taken  as  -vx  (Equation  l-*tU))  for  the  ij  element 
corresponding  to  the  beam  between  mass  1 and  3.  The  program  is  coded  so  that 
for  selected  ij  pairs  (specified  by  the  input),  the  internal  force  acts  only 
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on  mass  J,  not  on  mass  i.  In  Figure  1-lT,  beam  1-3  would  be  treated  in 
this  way,  with  its  internal  forces  acting  only  on  mass  3.  The  only  restriction 
is  in  the  numbering  of  the  masses;  the  mass  representing  the  torso  must  be 
assigned  a larger  number  than  the  mass  representing  the  seat. 

1.3.13  Integration  Routine 

The  equations  of  motion  derived  earlier  are  integrated  numericeLlly  in 

computer  program  KRASH.  The  integration  scheme  employed  is  a modified 

predictor-corrector  method.  Figure  I-I8  shows  a time  history  of  a typical 

response  quantity,  say,  u.  Assume  that  u and  u are  known  at  time  t and  all 

previous  times.  The  predicted  value  of  u at  t + At,  u is  computed  as 

^t+At 

u = u.  ..  + 2Atu  (1-126) 

Pt-At  ^ 

This  is  shown  in  the  Figure  I-I8  by  the  upper  dashed  line  of  slope 
Using  this  predicted  value  of  u at  t + At,  the  derivative  u at  t + is 


Figure  I-I8,  Numerical  Integration  Scheme 
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computed  from  the  equations  derived  in  the  theory.  This  derivative  is  then 
averaged  with  u at  t to  determine  a corrected  value  of  u at  t + At : 


u 

®t+At 


“t  * 2 


(1-127) 


This  is  shown  by  the  lower  dashed  line  in  Figure  l-l8.  The  final 
value  used  for  is  a weighted  average  of  the  predicted  and  corrected 

values.  For  the  present  program,  the  weighting  used  is  U to  1 in  favor  of 
the  corrected  value.  Hence 


= 0.8u 


+ 0.2u 


(1-128) 


A computational  flow  diagram  for  this  scheme  is  shown  in  Figure  1-19. 
Since  the  technique  is  not  self-starting  (due  to  the  u^  ^ term  in  equa- 
tion (1-126)),  a simple  Euler  integration  is  used  for  the  first  step.  The 
subroutine  DERIV  encompasses  all  the  equations  derived  in  the  theory,  and 
is  used  to  compute  the  derivatives  of  all  the  variables,  knowing  the  variables 
themselves . 

The  iteration  shown  by  the  dashed  line  in  Figure  1-10  is  often  employed 
with  this  scheme,  the  iteration  continuing  until  the  predicted  and  corrected 
values  agree  to  within  a specified  tolerance.  However,  the  iteration  is  not 
used  in  the  program  KRASH.  Due  to  the  large  number  of  computations  in  the 
subroutine  DERIV,  it  is  more  efficient  to  go  through  DERIV  only  once  per  time 
step,  and  to  choose  a sufficiently  small  time  step  At  to  obtain  the  desired 
accuracy. 

The  subroutine  IC  computes  all  the  initial  values  (at  t = 0)  of  the  vari- 
ables necessary  for  the  integration  to  proceed.  Section  1.3.15  discusses 
the  calculations  of  these  initial  conditions. 


I 
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Symmetrical  Modeling 

Program  KRASH  is  coded  to  allow  the  analysis  of  a half-airplane  model  for 
siTiunetrical  impact  conditions.  Figure  1-20  illustrates  a very  simple  model  of 
an  airplane,  in  which  masses  1,  2 and  3 represent  the  fuselage,  U suid  5 side 
fuselage  structure,  6 the  landing  gear  and  T the  wing.  If  the  airplane  itself 
is  symmetrical,  which  is  normally  the  case,  then  the  program  user  need  only 
input  data  for  masses  1 through  7.  Beam  data  is  input  for  the  following  ij 
pairs : 

1-2 
2-3 
2 - k 
U - 5 
U - 6 
5-0 
5-7 

The  notation  5-0  means  that  a beam  connects  mass  5 to  an  equal  mass  on 
the  opposite  side  of  the  airplane.  (Mass  9 is  not  input  to  the  model,  so  beam 
5-9  cannot  be  input.)  Program  KRASH  then  analyzes  the  response  of  the  model 
consisting  of  the  left  hand  side  of  the  airplane  only,  for  a symmetrical  impact 
condition.  A symmetrical  impact  is  one  containing  no  roll  or  yaw  suigles, 
no  lateral  eg  velocity,  and  no  roll  and  yaw  angular  velocities. 

For  the  half-airplane  analysis,  the  "missing"  half  of  the  airplane  is 
accounted  for  in  the  following  manner; 

• For  traverse  beams  such  as  5-0,  the  conditions  at  beam  end  J are 
set  equal  to  those  at  i,  with  appropriate  changes  in  sign,  x,  z euid 
6 motions  are  equal;  y,<j>and  ij*  motions  are  opposite. 

• The  X,  Z and  M internal  beam  loads  acting  on  a mass  on  the  plane  of 
symmetry  are  doubled  to  account  for  the  presence  of  opposite  side  beams 
that  are  not  explicitly  modeled.  Prior  to  this,  loads  for  beams  lying 
entirely  in  the  plane  of  symmetry  (such  as  1-2  and  2-3  in  Figure  1-20) 
are  halved.  This  is  done  so  that  when  the  X,  Z and  M internal  loads 
are  doubled  on  a mass  such  as  2,  the  contribution  from  beams  such  as 
1-2  will  not  be  doubled.  Y,  L and  N loads  are  set  to  zero. 
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The  input  data  defining  the  characteristics  of  beams  in  the  plane  of 
symmetry,  such  as  1-2,  should  include  the  contributions  from  both  the  left 
and  right  hand  sides  of  the  airplane.  The  same  is  true  for  mass  data  for 
masses  in  the  plane  of  symmetry,  such  as  1,  2 and  3-  In  other  words,  "half- 
airplane" data  is  not  utilized  for  beams  or  masses  in  the  plane  of  symmetry. 
While  this  convention  differs  from  that  used  in  certain  other  standard 
structural  analysis  computer  progreims,  it  is  consistent  with  the  f\ill  air- 
plane models  used  with  previous  versions  of  KRASH. 

For  unsymmetri cal  impact  conditions,  the  program  user  can  still  input 
a half-airplane  model.  For  this  situation,  the  program  internally  generates 
the  rest  of  the  model  and  performs  the  analysis  using  a fiHl-airplane  model. 
For  example,  the  program  would  internally  generate  the  data  for  masses  8, 

9,  10  and  11  in  Figure  1-20,  as  well  as  beams  2-8,  8-9,  8-10  and  9-11* 
Therefore,  the  program  user  can  always  work  with  a reduced  size  data  set 
representing  half  an  airplane,  as  long  as  the  airplane  itself  is  symmetrical. 
Any  change  made  to  the  half-airplane  data  set  will  automatically  be  included 
in  the  opposite  mass  or  beam  when  the  program  generates  the  full-airplane 
model. 


1.3.15  Initial  Conditions 

Initial  values  (at  t = 0)  are  required  for  the  following  quantities: 


*i*  ^i’  "i 


♦i.  0i. 


U^,  V^,  w. 


Pi, 


i = 1,2,. ..,N 


(1-129) 
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The  initial  conditions  which  are  input  to  the  program  are  the  following: 


X , y , i 
g’  'g  g 


P'.  q'.  r' 


vehicle  c.g.  velocity  In  ground  axes 

vehicle  angular  velocity  in  body  axes 

vehicle  attitude,  Euler  dhglefs  of  vehicle 
relative  to  ground 


First  the  total  vehicle  weight  is  computed  from  the  individual  mass 
weights  (input): 


W = 52  w 
TOT  ^ i 


(1-130) 


Next,  the  location  of  the  overall  vehicle  c.g.,  in  the  Hx"y"z"  coordi- 
nate system,  is  determined  as  follows: 


1 N 

x"  = — ^ 'y  w x" 


TOT  i=l 


1 N 

■ ' W S '^i^i 

^ TOT  i=l  ^ ^ 

1 N 

*0  = vT^ 

^ TOT  i=l  ^ ^ 


(1-131) 


An  expression  for  the  ground  coordinates  of  is  written: 


yA  + [A«]  /y' 


(1-132) 


r 


or 


va^ 


‘IJ  ■ ''*J  * •>  ■ 

where  x^,  y^,  and  are  currently  unknown,  and 


/ x'  \ / x" 

- 

I i 1 1 i 

G j 

1 y;  > = - 1 y'l 

f z'  1 1 z" 

\ / \ 

-'s) 

or 


vip^j  = -(vippj^j-vgpp^j)  J ® 1.2,3 


(1-133) 


I 
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The  vipp  are  the  input  coordinates  of  the  m. , and  the  vgpp  were  com- 
ij  1 ij 

puted  in  Equation  (1-131),  so  the  vip^j  are  known.  In  Equation  (1-132),  [A'J 
is  known  from  the  input  Euler  angles  4> ' , 6 ' , and  4*'  • The  only  unknowns  in 
Equation  (1-132)  are  x„,  y_,  z„  (defined  as  vg,).  Clearly,  the  y ground 

u u U J 

coordinate  of  the  initial  c.g.  position  is  arbitrary,  so  this  is  set  to  zero. 


It  remains  to  determine  the  initial  c.g,  position,  both  fore-aft  and 
vertically.  This  position  is  determined  within  the  program  so  that  the 
vehicle,  in  the  desired  attitude  with  respect  to  ground,  "fits  into"  the 
wedge  formed  by  the  ground  plane  and  slope  with  the  lowest  springs  .001  inch 
off  the  ground  and  slope.  This  situation  is  illustrated  in  Figure  1-21. 

Point  0 is  the  origin  of  the  ground  coordinate  system,  with  x positive  forward 
and  z positive  down.  Point  0 is  located  at  the  intersection  of  the  horizontal 
ground  plane  and  the  slope.  The  vehicle  eg  is  shown  in  Figure  1-21.  The 
initial  coordinates  of  the  c.g.  in  ground  axes,  shown  as  Xq  and  Zq  in  Figure  1-21, 
are  desired.  Also  shown  are  two  typical  masses  m^  and  mj , with  external 
springs  radiating  outward. 
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Vehicle  Initial  Position 


The  procedure  for  determining  and  z^,  is: 

1.  Write  equations  for  the  ground  and  slope  coordinates  of  the  end 
points  cf  all  the  external  springs  in  the  model,  assuming  the 
vehicle  eg  is  located  at  point  0. 

2.  Search  for  the  largest  magnitude  of  these  distances,  both  in 
ground  and  slope  axes. 

3.  From  this  information,  calculate  and  z^  to  give  .001  inch 
ground  clearance  at  the  spring  ends. 

li.  Finally,  calculate  x„. 

(j 

From  Equation  (1-69),  the  z coordinate  of  the  end  of  spring  ik  is 
known  as : 


(1-13U) 


From  Equation  (1-132)  va^^  is: 


vai3  = ^^^3  ^ 


1£ 


(1-135) 


i=l 


once  vg_(z_)  is  known.  Combining  (l-13^)  and  (1-135)  and  utilizing  (l-7l) 
3 O 


for  dvc 


ijk’ 


^"i3k  = ^«3  ^ I ^'31  ^ ^3K^: 

e=i 


ik 


(l-136a) 


The  X coordinates  are  obtained  similarly  as 


VC 


ilk 


= vg. 


f=i 


l£ 


+ A 


11 


ilk  ik 


(1-1 36b) 
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Nov  the  X and  z coordinates  of  all  the  spring  ends  are  calculated  using 
Equation  (1-136)  with  vg^  = vg^  = 0 (vehicle  c.g.  at  point  O).  The  zsl 
slope  eixis  coordinate  is  calculated  also  for  all  the  spring  ends.  This  is 
given  simply  by 


zs£  = sinp  (vc  ) + cosp  (vc  ) 


(1-137) 


The  results  of  equations  (l-136a)  and  (1-137)  are  now  searched  for  the 
maximum  magnitudes : 


largest  vc^^j^  for  all  ik  pairs 


!ot  = largest  zsl.,  for  all  ik  pairs 
max  ' 


(1-138) 


The  calculation  of  z^  and  Zg^  , shown  in  Figure  1-21,  is  now  possible. 
These  are  given  by 


Zf,  •=  -z  -.001 


(1-139) 


It  still  remains  to  calculate  Xg  from  Zg  and  Zg^^.  (zg  and  Zg^^  uniquely 

define  the  vehicle  c.g.  position,  but  we  want  x and  z , not  z„-  and  z.,  for 

vj  (j  SL  G 

use  in  equation  (1-132).) 

Obtain  first  the  projection  of  z„,,  on  the  vertical  Oz  axis  as 

oL 


-Z  /cosp 


(1-H*0) 


Then  Az  is  given  by 


' -"g  -"slg 


(1-lUl) 
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Finally  is  computed  as 


Xq  = dz/tanP  (l-ll*2) 

For  3=0,  don't  use  equations  (l-lUo)  through  (l-lU2),  but  set  x^  = 0. 

U 

For  3 = 90,  set  x = z-^  (both  are  negative).  Finally,  -modify  the 

O oJj 

calculated  values  of  x_  and  z_  by  the  input  values  x__„  and 

u vj  GIN  GIN 

- *GIN 
Equation  (I-IU2) 

(l-li*3) 

^G  ^G  ■ ^GIN 
Equation  (1-139) 

Positive  input  constants  x„^„  and  z„t.t  define  the  distances  aft 
and  up  in  ground  axes  that  the  initial  c.g.  position  is  relative  to  the 
"wedged"  position  shown  in  Figure  1-21.  Normally  these  are  input  as  zero, 
since  the  "wedged"  position  is  where  one  usually  wants  the  vehicle.  However, 
for  certain  special  cases  it  may  be  desired  to  drop  the  vehicle  a finite  dis- 
tance or  to  impact  entirely  on  the  slope.  It  should  be  noted  that  the 
z^  in  Equation  (I-IU3)  is  actixally  used  to  calculate  A z in  Equation  (l-li*l)  . 
This  means  that  if  a nonzero  Zqjjj  is  used,  the  vehicle  will  be  lifted  up  from 
the  "wedged"  position  and  then  moved  forward  until  almost  contacting  the 
slope.  Thus  one  can  start  any  distance  up  the  slope  to  get  clear  of  the 
ground  plane  if  so  desired.  The  slope  extends  indefinitely.  With  x^  and 
now  known  (and  y„  set  = 0),  we  can  use  Equation  (1-132)  to  calculate  the 
initial  values  for  all  the  mass  ground  coordinates. 

In  determining  the  initial  position  of  the  airplane,  only  the  mass  points 
having  external  springs  are  checked.  The  program  assumes  that  the  user  has 
assigned  springs  to  the  masses  at  the  lower  extremities  of  the  vehicle.  If 
thin  is  not  the  case,  the  program  may  place  the  vehicle  in  an  initial  position 
wherein  some  of  the  masses  are  below  the  ground  plane  (z  = 0).  This  will 
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r.ar.ifest  itself  by  negative  potential  energy  for  any  "subterranean"  masses. 

If  this  occurs,  the  user  should  correct  the  model  by  adding  external  springs 
to  the  lowest  masses. 

To  perform  the  preceding  calculations,  is  required,  which  is 

used  in  Equations  (l-136a,  b).  From  the  definitions  of  the  transformation 
matrices  (refer  to  Table  1-1  earlier): 

{Aj  - [A'][AJ]  (1-lUU) 

Lt_  1th  body  axes  to  c.g.  axes 
.g.  axes  to  ground  axes 


Equation  (l-l^ti*)  also  provides  a method  for  determining  the  initial  values  for 
*^1’  ®i’  '^i'  determined  from  ~ ~ sin0^  as 


6^  = -sin 


-1 


i31 


(I-1U5) 


Next  is  found  from  A.  = sini.  cos0.  , since  0.  is  now  known; 
1 i32  i i 1 


(I-IU6) 


Finally,  is  obtained  from  ~ sin'l'^: 


(1-1J*7) 


Note  that  equations  (l-l45)  throvigh  (I-IU7)  all  involve  arc  sine,  so  the 
difference  between  a positive  and  negative  angle  can  be  detected. 
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Next  the  linear  and  angular  velocities  are  calculated.  Differentiating 
equation  (1-132): 


(1-11*8) 


v'r.ere  zr.e  first  vector  is  input  and  the  second  vector  is  given  by 

iiu-’ticn  kl-133).  CA'D  is  treated  the  same  as  Ca.J  was  earlier  in  the  theory: 

U’]  = [a']  [d']  (1-11*9) 

where 


0 -r'  q’ 

r'  0 -P'  (1-150) 

-q'  p'  0 _ 

p' , q'  and  r'  are  input  constants  defining  the  initial  vehicle  angular 
velocity. 

Once  x^,  y^,  and  are  known  from  (I-1U8),  these  are  simply  transformed 
into  body  sixes  to  obtain  u^ , v^,  and  w^: 

(*i 

T ] . 

j^i 

("i 

The  angular  velocities  in  body  axes  (p^,  q^,  and  r^)  are  equal  to  the 
vehicle  angular  velocity  (p',  q',  and  r'),  transformed  from  c.g.  axes  to  ith 
mass  axes: 
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(1-151) 


Next  we  determine  6^ , and  “by  Equation  (1-115). 


where 


[Aj 


’l  sin()>^  tan6^  cos<|>j^  tanS^^ 


cos<(). 


-sin^. 


i 

1^0  sec6^  cos(|>^  secS^^ 


(1-152) 


(1-153) 


(1-15^) 


If  The  only  remaining  initial  conditions  to  be  determined  are  the  beam 

ij 

j(  deflections  and  forces.  At  present,  the  program  is  written  so  that  the  initial 

values  for  these  quantities  axe  zero.  This  is  only  true  for  the  case  where 
there  are  no  aerodynamic  loads,  so  that  all  masses  are  uniformly  accelerating 
downward  at  Ig.  Therefore,  if  aerodynamic  loads  are  used  in  the  program,  the 
initial  conditions  will  not  be  correct  since  internal  forces  would  be 
required  to  balance  the  aerodyneunic  and  gravity  forces.  However,  this  is  not 
considered  a serious  limitation  since  the  dynamic  impact  forces  in  a crash 
situation  will  overshadow  any  imbalance  present  due  to  nonzero  aerodynamic 
forces . 


Ir.  the  IC  suhrcutine,  the  overall  vehicle  moments  of  inertia  are 
calculated  from  the  individual  mass  inertia  data,  to  inform  the  user  as  to 
how  well  his  model  matches  the  overall  inertias  of  the  vehicle  being  analyzed. 
The  overall  inertias  are  given  by 


■ 2 ( 


I . + m . r . 1 
xi  i xii 


MM 

■ E * "i  V) 

i=l  . 

NM 

- E k,  * "i  4) 

* ^ ' 


(1-155) 


The  I I I , are  the  individual  mass  moments  of  inertia,  the  m,  are  the 
xi’  yi  zi  2 2 2 ^ 

individual  masses,  and  r . , r . , r . are  the  squares  of  the  distances  from 

xi  yi  zi 

the  eg  to  m^,  measured  normal  to  the  Hx"y"z"  axes  as  appropriate. 

4 - <yl  - yp^  * 


U-;  - * (z;  - 


U;  - xp^  * (y”  - yj)" 


(1-156) 
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1.3.16  Energy  Balance  Equations 

The  primary  objective  of  a crash  analysis  in  the  preliminary  design 
phase  is  to  determine  how  to  absorb  the  initial  vehicle  kinetic  energy  while 
maintaining  a livable  environment  for  the  occupants.  This  task  is  greatly 
facilitated  if  information  regarding  the  spatial  distribution  of  the  energy 
flow  through  the  vehicle  is  available.  With  this  objective  in  mind,  energy 
balance  equations  are  developed.  These  equations  do  not  alter  the  previously 
described  computational  procedures;  they  merely  provide  additional  information 
to  assist  in  understanding  how  the  initial  kinetic  energy  is  absorbed. 

The  total  system  energy  at  any  time  is  given  by  the  following  expression: 


^T  “ 


= KE  + PE  + SE  + DE  + CE  + FE 


(1-157) 


where 

Etot  “ Total  system  energy 

KE  = Total  kinetic  energy 

PE  = Total  potential  energy 

SE  = Total  strain  energy  absorbed 

DE  = Total  deunping  energy  dissipated 

CE  = Total  crushing  (externeil  spring)  energy  absorbed 

FE  = Total  friction  energy  dissipated  by  sliding  of 
external  springs  on  the  ground 

The  total  system  energy  remains  constant  during  the  analysis.  The 

total  kinetic  and  potential  energies  result  from  summing  the  energies  for 
each  mass  over  the  number  of  masses  (NM). 


KE 


NM 

E KE 
i = 1 ^ 


NM 

PE  * E PE. 
i = 1 ^ 


(1-158) 
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The  total  strain  and  damping  energies  are  obtained  by  siomming  the  strain 
and  damping  energy  for  each  internal  beam  element  (ij  pair)  over  the  NB  ij 
pairs . 

Nfi  MB 

SE  = S SE  DE  = E DE..  (1-159) 

IJ  = 1 iJ » 1 


from 

^rom 

over 


results  from  the  elastic-plastic  behavior  of  the  beam,  and  DE  results 

^ J 

its  damping  properties.  The  total  crushing  and  friction  energies  result 
summing  the  energies  for  all  the  individual  crash  springs  (ik  pairs) 
all  the  MSP  ik  pairs. 


FE  = 


(1-160) 


Referring  to  Equation  (1-15T),  at  time  zero  all  energies  are  zero  except 
KE,  PE  and  The  potential  energy  is  referenced  to  the  ground  plane. 

After  impact  with  the  ground,  the  kinetic  energy  decreases,  and  the  potential 
energy  may  increase  or  decrease.  Damping,  strain,  crushing  and  friction 
energies  all  increase  to  keep  Er|,Q^  constant. 

In  Equation  (1-159),  the  summations  over  the  iJ  pairs  exclude  those  iJ 
pairs  that  are  identified  in  the  input  as  DRI  elements.  This  is  done  because 
these  iJ  beams  and  their  masses  are  isolated  from  the  rest  of  the  system; 
the  forces  in  DRI  beam  iJ  drive  mass  J but  not  mass  i.  Also,  the  sximmations 
in  Equation  (1-158)  for  the  kinetic  and  potential  energies  exclude  mass  J 
in  a DRI  iJ  pair,  since  this  mass  is  isolated  from  the  system.  Thus,  if 
iJ  pairs  6-9  and  11-15  are  defined  as  DRI  elements,  the  summations  for  SE 
and  DE  will  exclude  these  iJ  pairs  and  the  summations  for  KE  and  PE  will 
exclude  masses  9 and  15- 

The  kinetic  energy  for  each  mass,  including  translational  and  rotational 
components,  is  simply 
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where 


I 4 1^X4 
xi  xyi  xzi 

I 4 X . 

xyi  yi  yzi 

I .X  , 1 . 
xzi  yzi  zi 


(1-162) 


(1-163) 


is  the  6x6  inertia  matrix  for  mass  i,  and  {vel^^}  is  a six-element 
vector  of  the  linear  and  angular  velocity  components  of  mass  i in  body-fixed 
axes. 

The  i^^  mass  potential  energy,  referenced  to  the  ground  plane,  is  given 


PEi  = -W^  z^ 


(I-I6U) 


Note  that  is  positive  downward,  measured  from  the  ground  plane. 

The  strain  energy  in  internal  beam  ij  is  computed  as  a continuous  summa- 
tion of  the  incremental  energy  contributions  from  each  integration  interval. 


previous 


(1-165) 
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F.„  dv. ^ is  the  internal  beam  force  (or  moment)  in  the  k direction,  mtilti- 
plied  by  the  incremental  beam  deflection  (or  rotation)  in  the  K ” direction  at 
beam  end  iM.  The  corresponding  J subscript  terms  apply  to  beam  end  JN.  The 
sum  of  these  constitute  the  incremental  strain  energy  for  the  integration  inter- 
val being  considered.  This  straightforward  formulation  automatically  accounts 
for  the  complexities  of  nonlinear,  coupled  deflections  and  unloading-reloading 

I 

behavior,  since  these  axe  considered  in  the  calculation  of 

The  damping  energy  dissipated  in  internal  beam  ij  is  computed  in  a sim- 
ilar meuiner: 


previous  I ■■ 


I th 

^Di?  internal  beam  damping  force  (or  moment)  in  the  i direction. 

The  crushing  and  friction  energies  are  calculated  as  a summation  of 
incremental  energy  changes.  These  are  first  expressed  in  terms  of  the  total 
external  spring  forces  acting  on  mass  m^,  and  the  incremental  displacements 
and  rotations  of  m^^.  These  total  forces  are  then  broken  down  into  the  com- 
ponents from  each  individual  external  sprir.j  (ik  pair),  and  further  into 
the  components  resulting  from  spring  loads  normal  to  the  surface  (crushing 
energy)  and  those  loads  in  the  plane  of  the  surface  (friction  energy). 

First  we  write 


jraj 


current 


(1-167) 


previous 


dqin^ 

drin. 

1 
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XCi*  NC^  are  the  six  forces  and  moments  acting  on  mass  i in  body- 

fixed  axes,  restating  from  all  ik  external  springs  attached  to  mass  i.  These 
are  given  by  Equation  (l-llO).  The  vector  jdx'^,  dy'j,  dz'j^,  dpin^,  dqin^^ , 
drinj^  } is  made  up  of  the  six  incremental  deflections  and  rotations  of  mass 
m^,  in  the  same  body- fixed  axes.  The  first  three  terms  of  this  vector  are 
given  by  a simple  rotation  transformation  of  the  incremental  deflections  in 
ground  axes,  which  are  obtained  directly  fron  numerical  integration  of  the 
equations  of  motion. 


dx' 

>'i^  = [*i] 

, dz’ . 


'dx. 
1 

>i 
,dz. 


(1-168) 


The  last  three  terms  of  the  incremental  displacement  vector  in  Equation 

(1-167)  are  the  incremental  changes  in  the  integrals  of  the  angular  velocities 

of  mass  m. , p. , q. , r . . These  are  the  incremental  rotations  of  mass  m.  in 
1 1 1 i 1 

body- fixed  axes  for  the  integration  interval  being  considered,  given  by  Equa- 
tion (1-12). 

The  negative  sign  in  Equation  ( 1-167 ) results  from  the  fact  that  the 
forces  acting  on  mass  m^  rather  than  the  forces  within  the  ik  springs  are 
being  considered.  A positive  deflection  of  spring  ik  resxats  in  a negative 
force  on  mass  m^.  The  energy  calculated  in  Equation  (l-l67)  includes  the 
energy  dissipated  by  the  sliding  of  spring  ik  on  the  ground  with  a friction 
coefficient,  since  the  forces  in  that  equation  include  the  ground  drag  loads 
due  to  friction. 

Now  separate  out  the  contribution  of  each  external  spring  (ik  pair),  by 
substituting  Equations  (1-109 ) and  (1-110)  for  the  [XC^,  YC^,  ...,  NC^J  in 
Equation  (l-l67).  This  reformulates  Equation  (I-I67)  into  a function  of  the 
individual  crash  spring  forces  FSP^^jj^,  where  i and  k refer  to  the  ik  spring  and 
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J refers  to  the  direction  of  the  forces  on  the  ik  spring.  These  forces  are 
shown  in  Figure  1-13.  The  resulting  equation  for  the  crushing/f riction 
energy  associated  with  each  spring  ik  is  the  following: 


current 


previous 


J = 1 


where 


I d vc^ 


(1-170) 


TERM.j^  = FSP.^^  (£il  + rx)  dqin.  - FSP^gl  ^ 

+ (rydrin^  - r^dqin^)  + FSP^^^rzdpin^  - FSP^^^rydpin^ 


k = 1 


= - FSP^22  ^^i2  dpin^  + F{5P^^2(^i2  * drin^ 

- ^2;dqin^  + FSP^^g  (rz.dpin^  - rxdrin^ ) + FSP^^^^xdqin^ 


k = 2 


"^^ik  “ ^^i23  dpin^  - (^13  + rz)  dqin^ 

k = 3 

^^il3  rxdrin^  + FSP^^^  (rxdqin^  - rydpin^ ) 


Finally,  Equation  (I-I69)  is  further  reduced  by  separating  the  individual 
spring  forces  FSP^^jj^  into  the  components  due  to  friction  and  crushing.  These 
are  denoted  by  and  FSPC^^j^,  respectively.  These  are  calculated  by 

using  Equations  (I-IO3)  and  {I-IO8),  with  XVOC^jj^  and  XVOCSL^jj^  broken  down 
further  into  crushing  and  friction  con5>onents.  This  results  in  the  following 
set  of  equations : 
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crushing 


friction 


crushing 


friction 


jjWOCC,.i=  ! 0 ( 


Contact 

with 

Horizontal 

Ground 


{ 


XVOCSLC 


ik 


i-|J, 


(1-172) 


XVOCSL 


i3k) 


. ( XVOCSL  ) 

|XV0CSLF^J  = XVOCSL. 


Contact 

with 

Inclined 

Slope 


crushing 

= - [Ai] 

T 

fxvocc.,  1 

1 ik  1 

} Contact 

= - [Ai] 

1 

. > with 

friction 

|fspFjJ 

1 

{XVOCF^^ 

l 1 Horizontal 

’ ) Ground 

crushing 

{fSPOjJ 

= - [Ai] 

T 

[abeta] 

|kVOCSLCj,^j  ] 

1 Contact 
f with 

friction 

{fSPFj^} 

11 

1 

T 

[abeta] 

jicvocsLr.^j  j 

[ Inclined 
I Slope 

(1-173) 


(CE,.  ) ^ = (CE,,  ) , Y'  FSPC 

ik  current  ik  previous  - ^ ^ 


ijk 


d VC, . + TERMC,, 
ij  ik 


(FE,J 


= (FE„ .) 


(1-17U) 


ik^  current  ik'  previous  Y FSPF  d vc  + TERMF 

j ^ ^ ij 


j = 1 


ik 


In  Equations  (I-I7U),  TERMC^j^  and  TERMF^j^  are  evaluated  using  Equation 
(1-171)  with  the  FSP's  replaced  by  FSPC  and  FSPF  as  appropriate.  Equations 
(1-17^)  give  the  crushing  and  friction  energies  for  each  spring  (ik  pair) 
used  in  equations  (I-I60)  to  calculate  the  total  CE  and  FE. 
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1.3.17  Stress  Equations 


The  program  is  coded  so  that  internal  heam  elements  can  be  analyzed  for 
stresses  and  checked  against  two  failure  criteria.  For  each  beam  element, 
the  location  of  the  maximum  bending  moments  in  two  planes  is  determined. 

This  occurs  at  one  end  of  the  beam  or  the  other,  never  in  the  middle.  The 
program  ceLLculates  the  stress  state  of  four  plane  elements  located  at  the 
four  sides  of  the  beam;  top,  bottom,  left  and  right.  Axial  stresses  due 
to  axial  load  and  bending  are  computed.  Shear  stresses  due  to  tension 
and  lateral  loads  are  also  calculated.  For  each  of  the  four  plane  elements 
of  a given  beam<^the  maximum  shear  stress  and  two  principal  stresses  are  then 
calculated.  These  are  used  to  form  ratios  of  actual  stress/failure  stress 
using  two  failure  criteria;  Maximum  Shear  Stress  Theory  and  the  Theory  of 
Constant  Energy  of  Distortion.  The  stress  calculations  in  no  way  alter  the 
scheme  for  confuting  internal  beam  loads  discussed  in  Section  1.3. 5. 3.  They 
are  an  additional  calculation  to  provide  the  program  user  with  a measure  of  the 
stress  state  for  the  beam  elements. 


Figure  1-22  illustrates  a typical  plane  element  as  viewed  from  the  top  (for 
the  top  and  bottom  elements)  or  from  the  left  side  (for  the  left  and  right  side 
elements ) . 


For  the  top  and  bottom  elements.  Figure  1-22  is  looking  in  the  direction  of 


the  positive  zb  beam  axis;  for  the  left  and  right  side  elements  in  the 

^ cl 


direction  of  the  positive  axis.  The  positive  sign  conventions  for  eucial 


stress  (T  and  shear  stress  a are  shown. 
X s 


The  axial  and  shear  stresses  are  calculated  as  follows : 


X = 


M z, 

y 1 

“T” 


Top 


(l-lT5a) 


s = 


QMx  + 


1.5  F 


X = 


M Z- 

_JL_i 
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1-5  F, 


Bottom 


(l-175b) 


«z  ^2 


X « 


= «Mx  _ 


z 

1.5  F 


Left  Side 


(l-175c) 


\ ^2 


Right  Side 


(l-175d) 


QM  1-5 

- X - z — 


and  are  the  sum  of  the  internal  loads  and  damping  loads 


F , F , F , M , 

X y z X 

given  by  Eqtiation  (1-58).  For  M and  M , the  maximum  value  at  either  point  i 

y z 


or  j is  used.  The  other  loads  are  constant  over  the  length  of  the  beam.  A, 


ly  and  are  the  average  beam  section  cross-sectional  area  and  moments  of 


inertia  about  the  ybj^j  and  zb^^j  axes,  respectively,  z^  and  z^  are  the  dis- 


tances from  the  neutral  axes  to  the  extreme  fibres  in  the  zb.  and  yb.  direc- 

i J 1 J 


tions,  respectively.  These  are  illustrated  in  the  following  sketch: 
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Q is  a sectional  shape  factor  relating  the  torsional  shear  stress  to 
the  applied  moment.  Reference  U,  Table  IX  contains  formulae  for  1/Q  for  a 
variety  of  commonly  encountered  structural  cross  sections.  The  1.5  factor 
for  the  direct  shear  contributions  in  Equations  (1-175)  results  from  the 
non-uniform  distribution  of  the  shear  stress  across  the  section.  The  peak 
value  at  the  center  is  used. 

After  the  axial  stress  it  and  shear  stress  a are  calculated  for  each 

X s 

c:’  the  four  locations,  the  maximum  shear  stress  <r,.  and  the  principal 

“max 

stresses  and  <t ^ are  calculated.  Reference  Table  II,  gives  the  fol- 
lowing equations : 


(l-lT6a) 

(l-176b) 
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These  stresses  are  used  to  calculate  "failure  ratios,"  which  are  the  ratios 
of  the  current  stress  to  the  failures  stress,  for  each  of  the  four  locations. 

The  failure  ratios  are  calculated  for  two  failure  theories;  Maximum  Shear 
Stress  Theory  and  Theory  of  Constant  Energy  of  Distortion.  For  the  former, 
the  maximum  shear  stress  from  equation  (I-I76)  is  divided  by  the  allowable 
shear  stress,  which  is  defined  internally  within  the  program  for  each  standard 
material  or  input  for  non-standard  materials. 

Failure  Ratio  = (1-177) 

shear 


For  the  Theory  of  Constant  Energy  of  Distortion,  Reference  U indicates  that 
elastic  failure  occurs  when  the  two  principal  stresses  satisfy  the  equation 


2(r 


2 

yield 


This  reduces  to 


2 ^2  2 
~ *^1  “'2  *^2  “ ‘^yleld 


(1-178) 


The  program  calculates  a reference  principal  stress  defined  as 

/ 2 F 

<^p  = ''  ""l  " ‘"1  ‘"2  ‘^2 

This  is  then  divided  by  'to  form  a Failure  Ratio  for  the  Theory  of 

Constant  Energy  Distortion: 

cr 

Failure  Ratio  = — E (1-179) 

‘Vleld 


Both  Failure  Ratios  from  equations  (1-177)  and  (1-179)  are  printed  out  at  each 
print  time  interval,  for  each  of  the  four  plane  element  locations  for  each 
internal  beam. 


SECTION  2 


PROGRAM  COirrROLS 


2.1  GENERAL 

This  section  contains  a brief  description  of 

• The  symbols  which  control  the  input-output  capabilities  of  KRASH, 

• the  program  subroutines , 

• a comprehensive  program  listing 

A User's  Guide  which  completely  defines  the  input-output  requirements  is 
provided  in  Volume  2,  Section  2. 

2.2  CONTROL  PARAMETERS 

Tables  2-1  through  2-h  list  the  various  input /output  control  parameters 
which  are  used  extensively  in  KRASH.  Table  2-1  presents  the  program  sizing 
constants,  and  shows  the  meocimum  allowable  sizes  for  the  different  types  of 
data.  Table  2-2  defines  the  program  execution  modes.  Tables  2-3  and  2-U 
define  the  print  output  and  plot  output  options,  respectively.  Table  2-5  and 
Figure  2-1  list  the  subroutines  used  in  KRASH  euid  provide  a brief  outline  of 
the  program  flow.  A Programmer's  Manual  (Reference  7)  contains  detailed  pro- 
gramming data. 


2-1 


J 


TABLE  2-1.  PROGRAM  SIZING  CONSTANTS 


CONSTANT 

MAXIMUM  VALUE 

DESCRIPTION 

NM 

80 

Nxjmber  of  masses 

NSP 

40 

Number  of  externeJ.  springs 

NB 

100 

Number  of  internal  beams 

NLB 

120 

Number  of  nonlineeu"  beam-direction  com- 
binations (KR  Tables) 

NHI 

8o 

Number  of  masses  having  non-zero  He^^^ , 

Hey^,  Hej,j^,  ^xyi*  ^ 

MVP 

— 

Reference  mass  number  for  volume  penetra- 
tion calculations 

NVCH 

5 

Number  of  volumes  for  occupiable  volume 
change  calculations 

NDRI 

100 

Number  of  DRI  beam  elements 

NMTL 

5 

Number  of  non-standard  beam  materials 

NACC 

50 

Number  of  input  acceleration  time-history 
tables 

NVBM 

100 

Niunber  of  internal  beams  having  non- 
standard maximum  deflections  for  beam  rup- 
ture. Standard  Value  = 100  (inches  of 
deflection  and  radians  of  rotation) 

NFBM 

100 

Nxmiber  of  internal  beams  having  non- 
standard maximvun  forces  for  beam  rupture. 
Standard  Value  = lElO 

NPH 

80 

Number  of  masses  having  non-zero  Euler 
angles  4>^",  G^"»  4^^" 

ND 

100 

Nvunber  of  internal  beams  having  damping 
ratios  different  from  that  specified  on 
card  0900 

NKM 

100 

' Number  of  internal  beams  for  which  the  full 
6x6  stiffness  matrix  is  directly  input 

NPIN 

100 

Number  of  internet!  beams  having  other  than 
fixed-fixed  end  conditions 

NNP 

50 

Niunber  of  node  points 

NUB 

100 

Number  of  unsymmetrical  beams 
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TABLE  2-2.  PROGRAM  EXECUTION  MODES 


1 

h 

I, 

t 

1 

I 


RUNMOD 

INPUT  DATA  SET 

DATA  SET  ANALYZED 

AIRPLANE  MODEL 

IMPACT  CONDITIONS 

0 

Full  Airplane 

Full  Airplane 

Unsymmetrical 

Unsymmetrical 

1 

Half  Airplane 

Half  Airplane 

Sy_  '.etrical 

Symmetrical 

2 

Half  Airplane 

F\ill  Airplane 

Synunetrical 

Unsymmetrical 

TABLE  2-3.  PRINT  OUTPUT  OPTIONS 


INPUT  CONSTANT 

CORRESPONDING  PRINT  DATA 

NSF 

Beam  strain  forces 

NTF 

Beam  total  forces  (strain  and  damping) 

ii 

NDE 

Beam  deflections 

NSPD 

Externeil  spring  loads  and  deflections 

'1 

NED 

Energy  breakdown  by  mass , beam  eind  spring 

l~] 

NS 

Beam  stresses 

1 

NRP 

Mass  displacement , velocity  and  acceleration  data 

1: 

i 

TABLE  2-5.  PROGRAM  KRASH  SUBROUTINES 


ERRATA  SHEET 

Report  No.  FAA-RD-77-189,  I 


PROGRAM  "KRASH"  THEORY 


1.  Delete  complete  Appendix  A,  Program  Listing,  and 
replace  with  enclosed  Appendix  A (Revised),  Program 
Listing,  pages  A“2  through  A“202. 

2.  The  enclosed  change  (yellow  sheet)  to  the  listing. 
Subroutine  'INPUT',  is  necessary  if  it  is  desired  to 
model  more  than  one  oleo  strut.  Refer  to  Page  A-130. 
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ERRATA  SHEET 

Report  No.  FAA-RD-77-189,  I 


PROGRAM  "KRASH"  THEORY 


lo  Delete  complete  Appendix  A,  Program  Listing,  and 
replace  with  enclosed  Appendix  A (Revised),  Program 
Listing,  pages  A“2  through  A"202. 

2.  The  enclosed  change  (yellow  sheet)  to  the  listing. 
Subroutine  'INPUT',  Is  necessary  If  It  is  desired  to 
model  more  than  one  oleo  strut.  Refer  to  Page  A“I30. 
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